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FOREWORD 


The work described in this report was performed by the United Technologies 
Research Center and by Pratt and Whitney Aircraft Group, Commercial Products 
Division, for the NASA-Lewis Research Center under Contract NAS3-22055. The program 
manager and principal investigator was Dr. Kevin P. Walker. The NASA-Lewis tech- 
nical project manager was Dr. Robert Kielb. 




Research and Development Program for Nonlinear 
Structural Modeling with Advanced 
Time-Temperature Dependent 
Constitutive Relationships 


1.0 SUMMARY 


The results of a twenty (20) month technical program, entitled '^Research and 
Development Program for Nonlinear Structural Modeling with Advanced Time-Temperature 
Constitutive Relationships”, are presented. This program was conducted by Pratt 
and Whitney Aircraft and the United Technologies Research Center for the NASA-Lewis 
Research Center under Contract NAS3-22055. The program included: (1) the evalua- 

tion of a number of viscoplastic constitutive models in the published literature; 

(2) incorporation of three of the most appropriate constitutive models into the 
MARC nonlinear finite element program; (3) calibration of the three constitutive 
models against experimental data using Hastelloy-X material; and (4) application of 
the most appropriate constitutive model to a three dimensional finite element analy- 
sis of a cylindrical combustor liner louver test specimen to establish the capa- 
bility of the viscoplastic model to predict component structural response. 

Based on a qualitative evaluation in Task I, three viscoplastic constitutive 
theories of material behavior, namely: (1) Walker’s functional theory; (2) Miller’s 

theory; and (3) Krieg, Swear engen and Rohde’s theory, were chosen for further study 
in Task IV. Evaluations of the three theories were based on a comparison of the 
predicted and the experimentally observed constitutive response of Hastelloy-X 
material developed in Task II. Each theory was incorporated into a subroutine 
of the MARC nonlinear finite element program in Task III. (A listing of the 
FORTRAN subroutine delivered to the NASA-Lewis Research Center is given at the end 
of the report.) Based on the evaluations in Task IV, Walker’s functional theory 
was chosen to analyze the structural response of a combustor liner induction rig ^ 
test specimen under temperature and strain cycles comparable to those in actual 
engine operation. The MARC subroutine incorporating the functional theory was 
delivered to the NASA-Lewis Research Center and installed for use on the UNIVAC 1100 
computer . 

The three viscoplastic constitutive theories were incorporated into the MARC 
program by means of an initial stress technique. All of the material nonlinearity 
in the constitutive equations is incorporated into an initial load vector and 
treated as a pseudo body force in the finite element equilibrium equations. 

The structural analysis of the combustor liner specimen under complex thermo- 
mechanical loading cycles has shown a deficiency in the currently proposed visco- 
plastic constitutive theories. These theories use various internal state variables 
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which employ a hardening-recovery format to determine their growth. The growth 
of the internal state variables due to hardening takes place only in the presence 
of inelastic deformation, while the recovery of the internal state variables can 
take place in the absence of inelastic deformation. Because internal state variable 
growth due to hardening can only take place in the presence of inelastic deformation, 
no change in the state variables (other than recovery or annealing) can take place 
during the elastic unloading phases of a thermomechanical cycle. Experimental evi- 
dence suggests that the state variables can change with temperature variations 
during the elastic unloading phase of a thermomechanical cycle. 
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2 . 0 INTRODUCTION 


f 


2.1 Program Need 

The overall operating cost of the modern gas turbine engine is greatly 
influenced by the durability of combustor and turbine structural components while 
operating at high temperatures. During each flight cycle, these components undergo 
large thermally induced stress and strain cycles which include significant amounts 
of creep and relaxation. Requirements for higher turbine inlet temperatures, uni- 
formity of temperature profiles, increased turbine cooling, and reduced emission 
levels have increased heat and pressure loads on combustor liners, while limiting 
the amount of air available for cooling the liners. Creep/fatigue cracking and 
creep buckling distortion of combustor liners reduces turbine durability through 
higher temperatures, and through impact and erosion damage caused by liberated 
pieces of the liner and hard carbon particles. 


2.2 Program Objectives 

The current program addresses a critical issue in the development of advanced 
life prediction technologies - the need to develop advanced viscoplastic constitu- 
tive relationships and the evaluate their effectiveness in predicting the cyclic 
stress-strain history at critical fatigue locations in major hot section components. 

In order to meet the objectives of the current program, the advanced consti- 
tutive models have been specifically calibrated against Hastelloy-X material which 
is used in the construction of combustor liners. But the constitutive models 
themselves, and the methods developed to evaluate the material constants required 
to support the models, are generally applicable to other materials. In addition, 
the modular form of the computer routines developed in the program should allow 
their incorporation, with necessary changes, into any nonlinear structural analysis 
finite element program. 


2.3 Program Relevance 

The current program provides an improved structural analysis capability and 
forms a sound base for future improvements in mission life analysis programs. With 
such screening tools available it will be possible to ensure that only promising 
new candidate schemes are put into hardware without waiting for one or two years of 
field operation to provide a feedback to the design scheme. 

Improvements in hot section component durability will have a profound impact on 
the operating costs of high bypass ratio turbofan engines through reduced maintenance 
cost. Equally important is the need to achieve reduced emissions levels and fuel 
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consumption while maintaining the durability of hot section components. The 
achievement of advancements in hot section component durability technology will ^ 

directly support the maintenance of U.S. leadership in the aircraft gas turbine 
engine industry. 


2 . 4 Background 

The overall operating cost of the modern gas turbine engine is significantly 
affected by the durability and efficiency of the major hot section components, the 
combustor and turbine. Primary responsibilities of the combustor are gas tempera- 
ture level and pattern control, required for efficient turbine operation, and ex- 
haust emission control at the various flight operating conditions. These goals 
are accomplished by the precise metering of air throughout the combustor structure. 
The high pressure and high combustion gas temperature characteristics of this 
environment require that the combustor liner be cooled for durability. 

These requirements for control of exit gas temperature, emissions, and metal 
temperature generate an intense competition for utilization of combustor airflow. 

The more aggressive performance, efficiency, and emission goals set for future en- 
gines emphasize the need for development of durable combustor structures which can 
operate with reduced levels of cooling air. This requires detailed knowledge of the 
operating environment and the ability to accurately predict structural response for 
these loadings. 

At high operating temperatures, the time dependence of the inelastic behavior 
of the structural material is a significant consideration in the design and analysis 
of combustor liners. Currently, most combustor liners are constructed of sheet 
metal louvers which rely on convective film cooling. Details of a louver geometry 
are shoyn in Fig. 1 where the individually formed pieces are seam welded together 
in the fabrication of the louver liner. Cooling of each individual louver is 
accomplished by compressor discharge air, introduced into the chamber created by 
the louver lip and knuckle region to generate an insulating film of cooling air on 
the downstream panel. Degradation of the cooling film, as it moves along the panel, 
results in increased metal temperature in the seam weld and louver lip regions. The 
thermal gradient between these regions and the knuckle of the next louver is a 
function of the engine operating condition. Figure 1 indicates the typical knuckle 
location (B) which is generally cooler than the typical louver lip location (A). 

The temperature difference between these two locations creates thermal stresses in 
the combustor liner. At the maximum power point, the thermally induced stress and 
strain state at the louver lip has significant amounts of plasticity which, when 
repeated over subsequent engine flight cycles, results in creep/fatigue failure of 
the liner. In large diameter annular combustor liners, fatigue cracks initiate in 
the louver lip region and grow axially toward the seam weld. Continued growth can 
result in distortion of the local cooling air stream, resulting in over— heating and 
burning of the liner and serious damage to the turbine components. Typical liner 
failures of this type are shoxm in Fig. 2. Analytical modeling and life prediction 
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of this kind of structural response has been hampered by lack of well-calibrated 
temperature data, high temperature cyclic material information, and precise engine 
failure data. With the advent of improved thermal and structural analysis capa- ? 

bilities and the development of test facilities to provide well-controlled component 
failure data, new structural modeling strategies and life prediction systems can be 
proposed and calibrated for the development of improved combustor liner structuires. 
With such analytical screening tools available, it will be possible to ensure that « 

only promising new candidate schemes are put into hardware without waiting for one 
or two years of field operation to provide feedback to the design scheme. 

t 

Life prediction of combustor liners requires a thorough knowledge of the thermal’ 
environment, accurate material characterization, general capability stress analysis 
techniques, calibrated failure data, and a life prediction model. The advent of 
general purpose finite element programs has provided a capability for detailed 
modeling of engine hot section components where the component structural response 
is dominated by significant amounts of nonlinear material behavior. However, the 
structural analysis of such components is usually carried out in nonlinear finite 
element programs, such as MARC (Ref. 1), which possess a degree of programming 
sophistication far in advance of their constitutive material modeling capability. 

Most nonlinear finite element programs, in keeping with an accepted practice 
in the elevated-temperature design community, partition nonlinear elevated-tempera- 
ture material behavior into rate-dependent "creep" and rate- independent "plasticity" 
components. Each component is assumed to obey a separate constitutive law in which 
the rate-dependent creep strain is assumed to be independent of the time- independent 
plastic strain and the plastic strain is assumed to be independent of the creep 
strain. No provision for interaction between creep and plastic behavior is present 
in the constitutive equations. These theories are therefore of questionable value 
at high temperatures where such interactions are known to occur. 

/In addition to the inability to model the interaction between creep and 
plasticity, most of the classical plasticity theories suffer from an inability to 
accurately model material behavior under cyclic load conditions. Under such 
loading conditions the classical theories are unable to predict the strain harden- 
ing/softening characteristics of the material. These limitations in constitutive 
modeling behavior have recently been discussed at a number of symposia (Refs. 2-6) 
and also in some depth by Krempl (Ref. 7). 

Sincef the constitutive material model lies at the heart of all general purpose 
nonlinear finite element codes, including the MARC program, it is clear that the 
development and efficient implementation of new viscoplastic constitutive models 
into such programs is needed before the full potential of such codes can be fully 
utilized for advanced high temperature component structural analysis. 
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State-of-the-art finite element methodology and life prediction algorithms 
have been evaluated against failure data from a combustor liner induction rig test 
specimen under contract NAS3-21836 (Ref, 8), with the NASA-Lewis Research Center. 
The work described in the present report provides an enhanced ability to assess the 
degree of cyclic nonlinear structural behavior generally associated with thermoj 
mechanical fatigue damage, A flow chart indicating where better material modeling 
capability fits into an improved life prediction scheme is shown in Fig. 3, 


2.5 Scope of Program ^ 

A twenty (20) month program was conducted to develop and evaluate a number of 
time-dependent constitutive theories suitable for the structural analysis of com- 
ponents which operate at elevated temperature. The work effort was organized into 
the following tasks. 

Task 1 - Constitutive Model Selection 

Ten viscoplastic constitutive models, available in the published literature, 
were examined to assess their potential capabilities and limitations in representing 
the thermomechanical response experienced by components (such as a combustor liner) 
which operate in a cyclic elevated temperature environment. Based on this examina- 
tion, three (3) theories were selected for further evaluation and incorporation 
into the MARC nonlinear finite element program. 

Task 2 - Materials Testing 

Specimen tests were conducted using a representative combustor liner material, 
Hastelloy-X, to determine the required temperature dependent material constants and 
to evaluate each of the three selected constitutive theories. These tests utilized 
smooth, uniaxial bar specimens of Hastelloy-X together with a closed loop servo- 
hydraulic INSTRON testing machine. The following tests were performed on the uni- 
axial bar specimens: 

(a) Fully reversed cyclic stress-strain tests for a range of strain rates and 
temperatures consistent with the operating conditions of a representative 
combustor; 

(b) Creep tests at various stress levels and temperatures starting from various 
points on steady-state hysteresis loops executed at constant strain rate 
under fully reversed strain controlled conditions; 

(c) Stress relaxation tests at various temperatures starting from zero initial 
strain values on steady state hysteresis loops executed at various strain 
rates under fully reversed strain controlled conditions. 
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Tests (a) , (b) , and (c) were used to determine the material constants for each of 
the constitutive theories at temperatures ranging from 427 °C to 982®C (800°F to 
1800'^F). 

A "faithful cycle" thermomechanical test having a prescribed strain-tempef a- 
ture history appropriate to that of a critical combustor liner fatigue location • 

(location A at the louver lip in Fig, 1) provided an experimental thermomechanical 
hysteresis loop for evaluation of the constitutive theories under complex loading 
cycles involving large changes in temperature and strain rate, ^ 

Task 3 - Incorporation of Theories into Finite Element Program 

Each of the three theories was incorporated into a FORTRAN user subroutine of 
the MARC nonlinear finite element program. An initial stress technique was used 
in which the incremental constitutive relationship used in evaluating the incremental 
initial load vector was determined by a sub increment method. An integral recursion 
relation was developed for Walker’s functional theory which has proved to be stable 
and accurate, even for large incremental time/ load steps. Miller’s theory and the 
theory of Krieg, Swearengen and Rohde have been formulated as integral theories, 
but have not, as yet, been incorporated in integral form into the MARC program. 

These theories have been integrated over each sub increment by means of an explicit 
Euler forward difference method. 

Task 4 - Comparison and Evaluation of Candidate Theories 

The three candidate constitutive theories were evaluated by comparing the pre- 
dicted uniaxial response of Hastelloy-X material with the experimental results 
generated in Task 2. Based on these evaluations the functional theory was selected 
to analyze the structural behavior of a simulated combustion chamber outerliner. 

The ^FORTRAN subroutines for both the differential and integral forms of the func- 
tional theory were then delivered to the NASA-Lewis Research Center and installed 
in the MARC program for use on the UNIVAC 1100 computer. Two demonstration problems 
were then executed to ensure correct functioning of the computer program at the 
NASA-Lewis Research Center. 
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3-0 TASK I - TIME-DEPENDENT CYCLIC 
PLASTICITY MODEL SELECTION 


3.1 Constitutive Theory Review 

Ten isotropic viscoplastic theories, listed in Appendices 1-12, were examined 
in order to assess their potential capabilities and limitations in representing the 
thermomechanical structural response experienced by a jet engine combustor liner. 

In order to obtain commonality in the theories, the equilibrium (rest or back) 
stress has been denoted by , the drag stress by K, and the inelastic strain by t 


In the isotropic theories selected for further evaluation, it is assumed that 
the inelastic response of the material is incompressible with cy^^ = 0. The formu- 
lation of the theories then implies that: (a) the response to a hydrostatic loading 

is linearly elastic, and (b) that the shear response of the material can be deter- 
mined from its uniaxial response. The inelastic incompressibility assumption can 
easily be relaxed if future testing shows that the multiaxial response cannot be 
determined from uniaxial behavior. From the limited amount of available experimental 
data, it would appear (Ref. 9) that the incompressibility assumption for the inelastic 
response is not satisfied. 

Some brief remarks on unified viscoplastic constitutive equations in general 
will now be given before proceeding to review each theory. 

All of the theories listed in the Appendices, with the exception of Valanis’ 
theory (Appendix 12) and the theory of Bodner, Partom and Stouffer (Appendix 7), 
allow the expression for the inelastic strain rate under monotonic loading conditions 
to be written in the form 


' /T r-l / cr-ii I 

C--6--f- = f (3.1) 

where f ^ is a function of the argument ((a-^^)/K) which differs amongst the theories 
and whose inverse is denoted by the function f, so that ff~^ = 1. Under rapid 
loading conditions the derivative terms in Eq. (3.1) become paramount and the 
response of the material aprpoximates to 6 ^ Ee. The viscoplastic formulations 
therefore exhibit an elastic response under infinitely fast loading conditions. 

An inelastic response under infinitely fast loading conditions could be incorporated 
into the theories by replacing Eq. (3.1) with 



(3.2) 
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in which a - !E is the instantaneous inelastic stress-strain response. A suitable 
growth law for 5 would then be required. 

Under monotonic loading conditions the tangent modulus of metals is much smaller 
than the elastic modulus, so that for large strain values a/E << e, and the asymptotic 
stress-strain relation in Eq. (3.1) may be written as 


(3.3) 


^ +Kf(^)| as € — ± 00 . 


Cyclic hardening of the hysteresis loops can be accommodated by allowing the ^ 

equilibrium stress Q and the drag stress K to increase with an increase in the cumu- 
lative inelastic strain incurred by the material under cyclic loading conditions. 

Two repositories exist for incorporating cyclic hardening into the viscoplastic 
theories. First, consider the case where the equilibrium stress Q rapidly saturates 
to a limiting value at large strains, so that Q has reached a constant value at the 
tip of the first cyclic hysteresis loop reversal. Cyclic hardening can be incor- 
porated into the theories by assuming that the drag stress K increases with cumula- 
tive inelastic strain under cyclic loading conditions. If the material exhibits a 
positive strain rate sensitivity the cyclic stress amplitude increases with strain 
rate under strain-controlled conditions and the term Kf(e) increases with increasing 
values of e. At elevated temperatures and high strain rates the term Kf(e) in 
Eq. (3.3) may be the dominant term at the tips of the hysteresis loops. Under such 
conditions the hysteresis peak stresses may be approximated by the expression 
o ^ ± If the drag stress increases from an initial value of to a final 

saturated value of under sustained cycling, the fractional increase in stress 
amplitude, viz. | | / where Op = Kpf(e) and = Kjf(e), is independent of 

strain rate. The same relative amount of cyclic hardening is observed at all strain 
rates when the strain rate is large. Hardening of the drag stress K which is a 
scalar quantity, with cumulative deformation, introduces isotropic hardening into 
the viscoplastic formulation. At high strain rates the original stress amplitude 
±{oj ~ Kjf(e)} increases equally in both the compressive and tensile directions # 

to a final saturated value of ±{ap = Kpf(e)}. An example of this type of isotropic 
hardening is shown in Fig. 4. 

Another type of cyclic hardening can be incorporated into the viscoplastic 
theories by assuming that the drag stress rapidly saturates to a limiting value at 
large strains (or remains constant) . Cyclic hardening is then assumed to occur due 
to increases in the value of the equilibrium stress tensor at the peaks of the 
hysteresis loops due to sustained cyclic loading. At the hysteresis loop peaks the 
term Kf(e) is now constant for a given rate of straining and the stress amplitude 
increases due to increases in the magnitude of Q, This type of kinematic hardening 
produces a response quite different to that produced with an isotropic hardening 
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mechanism. Since the hardening occurs in the equilibrium stress Q, the relative 
amount of observed cyclic hardening will decrease with increasing strain rate. At 
sufficiently high strain rates where the term Kf(e) in Eq. (3.3) is paramount, 
virtually no cyclic hardening in the stress amplitude will be apparent. Kinematic 
hardening of the equilibrium stress produces the type of hardening observed in 
Fig. 5. In this figure, the peak stress amplitude increases with cycle number 
while the "yield stress" decreases with cycle number. This results in a plastic 
tangent modulus which increases with cycle number. An explanation of this effect 
can be obtained from an examination of Eq. (3.1). The production of an observable 
value of the inelastic strain rate, c, depends on the magnitude of the stress 
difference As the magnitude of ^ increases at the hysteresis peaks due to 

cyclic hardening, the magnitude of stress a required to produce "yielding", or 
an observable value of c in the reverse direction, decreases. In general, it is 
appropriate to incorporate both forms of cyclic hardening into the constitutive 
formulation. The hardening/softening mechanisms adopted by the unified viscoplastic 
theories are displayed in Table 1. 

In addition to differences in the functional form adopted for the grovjth of 
the inelastic strain rate with the state variables Q and K in Eq. (3.1), the visco- 
plastic theories also differ in the functional form assumed for the growth laws 
which determine the evolution of the state variables Q and K. These state variables 
are assumed to evolve according to differential constitutive equations which exhibit 
a hardening /recovery format. At large strain values these state variables saturate 
to limiting functional forms which, in the theories of Walker, Chaboche and Lee & 
Zaverl, are independent of strain rate at high rates of strain. In the remaining 
state variable theories the limiting functional forms of Q and K depend on the 
strain rate even at high rates of strain. 

Saturation of the equilibrium stress and drag stress state variables is 
achieved when the antagonistic effects of hardening and recovery cancel each other 
in the state variable evolution equations. If these equations are written in the 
form 




(3.4) 


and 


K = f4lcl-f5Klcl-fgK , 


(3.5) 


where fj^, £29 •••^5 functions of Q and K, and c is the inelastic strain rate, 
then the limiting saturation values are achieved when the growth rates of 0, and K 
become small. One may therefore set ^ = 0 and = 0 to obtain the limiting forms 
of 9, and K as 
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TABLE 1 


4 

MECHANISM FOR CYCLIC HARDENING/ SOFTENING IN THE 
UNIFIED VISCOPLASTIC THEORIES 


Hardening/ Softening^ 
due to both Q and K 

Hardening / So f t ening 
due to K** 

Hardening/ Softening 

^ ^•k-k'k 

due to 0, 

No Hardening/ 
softening 

Walker 
Chaboche 
Lee & Zaverl 

Miller 

Cernocky and Krempl 

Krieg, 

Swear eng en & 
Rohde 


Hardening/softening of the Q and K state variables implies a combination of 
cycle dependent kinematic hardening with an isotropic hardening component 

Hardening/softening of the K state variable implies a combination of ideal 
kinematic hardening with an isotropic hardening component 

Hardening/softening of the state variable implies cycle dependent kinematic 
hardening with no isotropic component 

No hardening/softening of the state variables implies ideal kinematic hardening 
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and 


a-- 


f4lcl 

K = 


(3.6) 


fslcl + fe ■ j 

These limiting forms depend on the strain rate at low rates of strain. 

At high rates of strain f 2 l^| aiid f^|c| >> f^. The state variables , 

then achieve the rate independent forms 


il = ±{f,(fi.K)/f2(n,K)} 


(3.7) 


and 


K={f4(n.K)/f5(n,K)} 


(3.8) 


The terms f 2 ^^|cl and f 5 Klcl in Eqs. (3.4) and (3.5) which are responsible for these 
rate independent limits may be called dynamic recovery terms, since they are acti- 
vated only in the presence of an inelastic strain rate. The terms f^l^^l and f^K 
in Eqs. (3.4) and (3.5) represent the effects of static thermal recovery and are 
activated even in the absence of an inelastic strain rate. If the dynamic recovery 
terms are absent the limiting values of the state variables may be written in the 
rate dependent forms 


c{f,(fi,K)/f3(n,K)} (3.9) 

and 

K = lcl{f4(fi,K)/f6(n,K)}- (3.10) 


The limiting state variable behavior of some of the viscoplastic formulations is 
listed in Table 2. It should be emphasized that although some formulations exhibit 
a rate independent state variable behavior at high strain rates, the viscoplastic 
theories still exhibit a rate dependent stress-strain response through the presence 
of the function f(e) in Eq. (3.3). 
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The behavior of the viscoplastic theories under creep and relaxation conditions 
is governed by the growth of the state variables. Figure 6 shows the steady state 
stress-strain hysteresis loop GJAHCFIG executed at a constant strain rate magnitude 
under fully reversed strain controlled conditions, while the equilibrium curve is 
shown as the "S” shaped hysteresis loop DEF. If, under steady state cycling, the 
stress is held constant at point A on the hysteresis loop the material will creep 
along the line ABC at constant stress. Initially, when the stress is held constant 
at point A, the equilibrium stress is determined by the stress value at point D. 

As the material creeps along ABC the equilibrium stress grows along the curve DEF 
according to the growth laws for Q in the Appendices. The difference between the 
actual stress and the equilibrium stress decreases to a constant value during creep. 
Since the creep rate depends on this difference, the creep rate also decreases 
(primary creep) and then remains constant (secondary creep). If the initial stress 
is large enough, so that the initial value of the equilibrium stress is at its sat- 
uration value, then the difference o-Q will remain constant in the creep test and 
secondary creep will occur without a primary transient. In the case where creep 
initiates at point J in Fig. 6, where the stress lies below the peak of the equili- 
brium stress-strain curve, no secondary creep takes place and creep terminates at 
point D on the equilibrium stress-strain curve. 

If the stress is held constant at point C on the unloading branch of the 
hysteresis loop, the initial value of a-Q ^ CF is small, and very little creep 
is observed initially. At elevated temperature the equilibrium stress Q will 
gradually decrease due to static thermal recovery (if included in the constitutive 
formulation), so that the creep rate gradually increases due to the increase in 
o-Q* Eventually, the decrease in the equilibrium stress is halted when the harden- 
ing of Qy governed by the fj_ term in Eq. (3.4), balances the recovery terms con- 
taining f 2 and f^. Secondary creep then takes place at the point C. The creep 
behavior on the unloading branch of the hysteresis loop may be contrasted with the 
creep^ behavior, at the same stress level, on the loading branch of the hysteresis 
loop. Here, the initial value of o-l^ is large, and a rapid primary transient in 
the creep curve is observed. This difference in creep response on the loading 
and unloading branches of the hysteresis loop is shown in Fig. 7. 

The behavior of the constitutive theories in relaxation is analogous to the 
creep behavior. Under stress relaxation conditions at constant strain the relaxa- 
tion rate depends on the difference o-Q so that if stress relaxation commences at 
point A in Fig. 6, the amount of relaxation is large. However, if stress relaxa- 
tion commences at point C on the unloading branch, very little relaxation is 
initially observed. As with the creep behavior, the hesitation in the relaxation 
curve observed on the unloading branch of the hysteresis loop eventually gives 
way to more rapid relaxation as the term o-Q increases due to thermal recovery 
(if included) of the equilibrium stress. 
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When the initial point at which creep or relaxation commences on the unloading 
branch of the hysteresis loop is below the equilibrium curve, at point I in Fig, 6, 
the theories predict that the stress will initially relax upwards towards the 
equilibrium curve during a stress relaxation test in which the total strain is held 
constant, or will creep in the negative direction (contract) at the positive stress 
value at point I during a creep test. The negative creep and negative (upward) 
relaxation will terminate on the equlibrium curve if the static thermal recovery 
containing f-^ is absent in Eq. (3.4) governing the growth of Q. However, if static 
thermal recovery is included, the equilibrium stress 0 . will eventually drop below 
the actual stress a at point I during a creep test and, with a-Q now being positive, 
the material will commence creeping in the positive direction. Similarly, the 
upward relaxation will cease when the equilibrium stress drops below the actual 
rising stress value and normal (downward) relaxation then ensues. Figure 8 shows 
that negative stress relaxation commences in compression at -21 MPa (-3 ksi) on 
the unloading hysteresis branch and progresses through zero to tensile values as 
the stress relaxes towards the equilibrium curve. 

With the foregoing remarks in mind a review is now given of ten viscoplastic 
constitutive theories available in the published literature. 


3.2 Walker’s Functional Theory 

This theory (Ref. 11) was developed in an integral form by modifying the 
constitutive relation for" a three parameter viscoelastic solid. Both the integral 
and the differential forms of the theory are summarized in Appendices 1 and 2. Two 
state variables, and K, are introduced into the viscoelastic theory to account 

for the effects of viscoplasticity. The equilibrium (rest or back) stress 
introduces nonlinear kinematic hardening into the model to account for the Bauschinger 
effect, while the drag stress K introduces isotropic hardening into the model to 
account for cyclic hardening or softening of the material. 

The growth law for the equilibrium stress contains both dynamic recovery and 
static thermal recovery terms. At high strain rates the thermal recovery term 
becomes insignificant in comparison with the dynamic recovery term and the 
equilibrium stress becomes independent of strain rate. In the growth law for the 
drag stress, static thermal recovery terms have been omitted, which permits the 
drag stress to be written in the integrated form shown in Eq. (3) of Appendix 1. 

This form has been found adequate in the modeling of Hastelloy-X behavior, but 
future applications may require the inclusion of static thermal recovery in the 
drag stress evolution law. 

The theory is capable of modeling the cyclic hardening and softening of 
hysteresis loops without the use of a yield surface. Material constants required 
to model cyclic hardening/softening are obtained from cyclic hysteresis tests so 
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that cyclic hardening and softening can be accurately modeled. Both the 
equilibrium stress Q and the drag stress K contribute to the cyclic hardening 
in the theoretical formulation. 

Creep, relaxation and strain rate effects are modeled by a power law (Eq. (^1) 
of Appendix 2) for the inelastic strain rate. This expression has been found 
adequate for the representation of strain rate effects encountered in a combustor 
liner material under service conditions where strain rates may vary from 10 ^ sec ^ 
to 10”^ sec“^. However, it appears that modifications are necessary if strain 
rates greater than about 10”^ sec*“^ are encountered. In particular, at higher 
strain rates, the power law expression for the inelastic strain rate predicts values 
of stress in a constant strain rate tensile test which are too large (cf. (Refs. 

10 and 20)). 

The integral expression for the equilibrium stress in Eq. (2) of Appendix 1 
contains the material constant n 2 which appears outside the integral. On differen- 
tiation with respect to time this integral expression reduces to the differential 
form given in Eq. (2) of Appendix 2. The terms containing 0 allow the equilibrium 
stress to change during nonisothermal elastic excursions in which the inelastic 
strain rate is zero. Reasons for the inclusion of the © terms in the growth law 
for the equilibrium stress are discussed in Section 6 of this report. 

Calculations may be made with the theory to construct a theoretical yield 
surface. For a given point in stress space one may construct a surface surrounding 
this point such that in going from the given point to each point on the surface 
the cumulative inelastic strain is the same for each point on the surface. If 
this theoretical construct is called the yield surface, then the theory allows the 
yield surface to translate in stress space by virtue of the presence of the equili- 
brium stress (kinematic hardening) and to expand in stress space due to the presence 
of the drag stress (isotropic hardening) . No provision for a rotation of the yield 
surface in stress space exists in the current theoretical formulation. All of the 
material constants in, this theory are functions of temperature and must be experi- 
mentally determined at each temperature of interest. 


3.3 Chaboche’s Theory 

This theory was developed by Chaboche (Refs. 12, 13), and is summarized in 
Appendix 3. The theory is similar in form to the preceding functional theory 
except that a yield surface concept is used. Inside the yield surface it is assumed 
that no inelastic deformation can take place. The use of a yield surface permits 
isotropic hardening to be modeled by an increase in the size of the yield surface 
rather than by an increase in the drag stress state variable, K. Hence, in this 
theory, K is assumed to be constant and the one-dimensional form of the inelastic 
strain rate may be written in the form 
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Icr-fll-Y-k 

K 


(cr-n) 

Icr-fil 



(3.11) 


where the brackets < > are defined by the relation 


<a> 


0 if a ^ 0 
a if a>0 . 


(3.12) 


Here, the isotropic state variable, Y, together with the yield surface concept, 
replaces the drag stress, K. Initially, the variable Y is assumed to be zero and 
inelastic deformation occurs only when |a-f^|>k. As Y grows with inelastic deforma- 
tion, the yield surface expands and inelastic deformation takes place only when 
la-^|>k+Y. At high rates of uniaxial strain and large strain values (where a 0) , 
Eq. (3.11) takes the approximate form 


a«n+Y+k + K€ (e— + CO)- 


(3.13) 


At high strain rates Q and Y become independent of strain rate, so that for very 
high strain rates 


CT 


l/n 


(3.14) 


This equation shows that at very high strain rates, no cyclic hardening is observed. 
Significant hardening will be observed only at lower strain rates when (^^ + Y + k) 
becomes comparable in magnitude with K e l/i^. 

The growth law for the equilibrium stress contains a 0 term which allows the 
equilibrium stress to change with temperature during nonisothermal elastic excur- 
sions. Chaboche’s theory. Walker’s theory, and Cernocky & Krempl’s theory, are 
the only viscoplastic formulations which permit such changes in the state variables 
to occur. 

The yield surface can translate and expand in this formulation, but it cannot 
rotate. All of the material constants are functions of temperature and must be 
experimentally determined at each temperature of interest. 
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3.4 Miller’s Theory 


This theory (Ref. 14) is summarized in Appendices 4 and 5. Instead of a power 
law for the inelastic strain rate, the theory uses a hyperbolic sine function. In 
addition, only the material constants Kq and 0’ are temperature dependent, Kq being 
the initial value of the drag stress K. 

For mono tonic one-dimensional loading. Miller’s theory may be written in the 

form 


C = 6 - 






(3.15) 


n=H,c-H,Be' {sinh(A,ri)}", 


(3.16) 


K=H2lC|(C2+ii A2K^/ A|)- H2C2B0'|sinh(A2K^)|"' 


(3.17) 


The equilibrium stress is assumed to harden linearly with inelastic strain. 

Only static thermal recovery is included in the growth law for the equilibrium 
stress. In the absence of a dynamic recovery term, the equilibrium stress continues 
to increase with increasing strain rate; at high strain rates the growth law does 
not become independent of strain rate. In modeling the stress-strain behavior of 
materials, this theory exhibits the characteristic tri-linear curve in Fig. 9. The 
first portion of the tri-linear curve in Fig. 9 corresponds to an "elastic" loading 
phase in which the inelastic strain rate and rate of growth of the equilibrium stress 
is small. On entering the inelastic region the equilibrium stress grows linearly 
with inelastic strain which gives rise to the second portion of the tri-linear stress- 
strain curve. This may be seen by noting that in this region the growth law for Q 
(with the neglect of the static thermal recovery term) may be written as 
0. ^ Hj^c , or since a << Ee, as « H^e. In this region, where the static thermal 
recovery term is small, the equilibrium stress increases linearly with strain. 

Since a << Ec, Eq. (3.15) governing the growth of inelastic strain, may be written 
for monotonic loading in the form 


0*^52 H + K 



(3.18) 


Usually, the drag stress K grows slowly with deformation and the rate dependent 
term in Eq. (3.18) may be assumed to be constant (for i = constant). The stress 
then increases according to the relation a + constant. Finally, the equili- 

brium stress reaches values sufficient to activate the static thermal recovery 
term in Eq. (3.16). This term increases very rapidly according to a hyperbolic 
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sine law and when it reaches a value equal in magnitude to that of the linear 
hardening term, the equilibrium stress reaches a constant saturated value. This 
gives rise to the third tri-linear portion of the stress-strain curve in Fig. 9. 

The equilibrium stress saturates rapidly with deformation in this theory so ^ 
that no provision exists to model cyclic hardening or softening by means of cyclic 
changes in the equilibrium stress. Cyclic changes in the stress-strain response 
are modeled by cyclic changes in the drag stress. In the evolution equation for 
the growth of the drag stress both dynamic recovery terms and static thermal 
recovery terms are present. At fast strain rates only the dynamic terms (those 
multiplied by |cl) are significant. Now at large strain values the saturated value 
of the equilibrium stress U may be determined from Eq. (3.15) by setting ^ = 0. 
Hence, denoting the saturated value with the subscript s, the saturated equilibrium 
stress may be written as 


iis* A| sinh I , (3.19) 

from which it is evident that the saturated equilibrium stress increases with strain 
rate. From Eq. (3.17), only the first bracket is paramount at fast strain rates. 
Within this bracket both Q and K increase with strain rate. When the drag stress , 
reaches a saturation limit at a given constant strain rate, one may set K = 0, and 
by neglecting the constant C 2 fri comparison with Q. and K, one obtains at satura- 
tion 



( 3 . 20 ) 


This relation was derived under the assumption of rapid strain rate. But from the 
construction of Eqs. (3 .15) - (3 . 17) , the relation may be shown to be valid at all 
strain rates. For example, if Eq. (3.20) holds at all strain rates for the saturated 
values of Q and K, then setting Eq. (3.17) equal to zero at lower strain rates 
where static thermal recovery is important produces the relation 


lcUB0'|sinh(A2Kf (3.21) 

Moreover, by setting Eq. (3.16) equal to zero, the saturated equilibrium stress is 
governed by the relation 


C=B0'-[sinh (A|iis)}"- 


( 3 . 22 ) 
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The relation between the saturated equilibrium and drag stresses given in 
Eq, (3.20) is now seen to hold at all strain rates by equating equations (3.21) and 
(3.22). One further relation between the equilibrium stress, drag stress and actual 
stress may be obtained by equating the saturated inelastic strain rates given in 
Eqs. (3.21) and (3.22) with that given in Eq. (3.15). One obtains ^ 




(3.23) 


The assumption is now made that under steady-state creep conditions, the inelastic 
strain rate in Eq. (3.15) is governed by the Garofalo equation, viz.. 


C=B0 


|sinh( Ao-g)} ■ 


(3.24) 


Hence, under steady state conditions 



(3.25) 


The saturated equilibrium stress is linearly related to the saturated stress while 
the saturated back stress depends on the one- third power of the stress. This one- 
third power relation for the saturated drag stress, which is derived from the 
effect of warm working at various temperatures on the subsequent room-temperature 
yield strength of aluminum and 304 stainless steel, was assumed by Miller in 
setting up the framework of the viscoplastic theory. 

Since cyclic hardening is incorporated into the isotropic drag stress state 
variable, hardening may be observed at all strain rates. At high strain rates 
the relative amount of cyclic hardening is the same for all strain rates as 
emphasized in Section 3.1. 

The state variables Q and K do not contain 0 terms so that during nonisothermal 
"elastic" excursions the state variables "stick" at the values they had on entering 
the "elastic" region. If a theoretical yield surface is constructed with this 
theory, the yield surface can expand and translate in stress space, but no rota- 
tion of the surface is permitted. 

The material constants n, Hj^, H 2 , A 2 , C 2 , and K of Appendix 4 are 

assumed to be independent of temperature. Those which depend on temperature are 
the elastic constants X and y; the initial value of the drag stress K^, which 
depends on the initial temperature at time t = 0; and the constant 0’ defined in 
Eqs. (7) and (8) of Appendix 4. 
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3.5 Lee & Zaverl’s Theory 


This theory (Ref. 15) employs a yield surface, which, in its general anisotropic 
form, can expand, translate and rotate. In this respect, it is the only theory^ 
considered in this report which allows the yield surface to rotate. However, 
in the more restricted isotropic form of the theory given in Appendix 6, no yield 
surface rotation due to deformation can occur. 

Cyclic hardening has been incorporated into the yield stress state variable Y, 
the drag stress state variable K, and the equilibrium stress state variable 
Since the inelastic strain rate in Eq. (1) of Appendix 6 depends on these state 
variables, it is appropriate to call them primary state variables. These state 
variables reach a saturation value at large strains by virtue of dynamic recovery 
terms alone, since no static thermal recovery terms are included in the theory. 

The dynamic recovery terms, which allow the primary state variables to reach their 
saturation values, contain the state variables a?j and Y^. These additional state 
variables may be called secondary state variables since they govern the growth of 
the primary state variables. Because static thermal recovery terms are absent in 
the growth laws for the primary state variables, the saturated values of Y, K and 

are independent of strain rate. Moreover, the omission of 0 terms in the state 
variable growth laws prohibits changes in the state variables during nonisothermal 
"elastic'* excursions. The theory is more demanding on computer storage requirements 
than the preceding theories due to the presence of the secondary state variables. 


3.6 Bodner, Par tom and Stouffer's Theory 

This formulation (Refs. 16-18), which is summarized in Appendix 7, differs 
from the other visoplastic theories by rejecting the concept of the equilibrium 
stress state variable It ts not, therefore, possible to describe negative 

creep effects at positive stress values, negative relaxation (see the experimental 
data in Fig, 8), or an aftereffect type of recovery in the material when the stress 
is removed. 

The theory is capable of modeling kinematic and isotropic hardening, but 
accurate modeling of cyclic effects will probably require the constant q (which is 
equal to one for purely isotropic hardening and equal to zero for purely kinematic 
hardening) in Appendix 7 to be taken as a function of the plastic work. 

Bodner and Stouffer have formulated a general anisotropic theory in Ref. (17), 
but only the isotropic form of this theory is shown in Appendix 7. Further comments 
concerning a generalization of the anisotropic model by Young may be found in 
Ref. (18). 

The theory does not contain 0 terms in the evolution equation for the state 
variable z. Consequently, this state variable cannot change during nonisothermal 
"elastic" excursions which precludes an accurate analysis of thermomechanical 
cycling. 
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If a theoretical yield surface is constructed from the model the surface can 
expand and translate in stress space but no rotation is permitted in an isotropic 
material. 


3.7 Krieg, Swearengen and Rohde's Theory 

This theory uses a power law to model the inelastic strain rate and was 
developed in a differential form (Ref. 19) by the authors. Both differential and 
integral forms of the theory are presented in Appendices 8 and 9. 

An equilibrium stress and back stress are included in the formulation but no 
explicit provision exists to model cyclic hardening in either of the state variables. 
That is to say, the constants A^ - A^ appearing in the growth laws for the equili- 
brium and drag stresses do not explicitly depend on the cumulative inelastic defor- 
mation. Cyclic hardening could be modeled by assuming that the constants A^ and 
A^ are small in the growth law for the drag stress, so that the drag stress grows 
slowly with deformation. 

The state variables saturate at large strain values when the static thermal 
recovery terms balance the linear strain hardening terms. No dynamic recovery terms 
are present in the formulation so that the saturated values of the state variables 
depend on the strain rate. Moreover, the linear hardening terms and the static 
thermal recovery terms in the state variable evolution equations produce stress- 
strain curves and hysteresis loops which exhibit the same tri-linear character as 
Miller's theory. The rapid final flattening of the stress-strain curve occurs 
when the static thermal recovery term grows rapidly with equilibrium stress accord- 
ing to the exponential law in Eq. (2) of Appendix 8. The sharpness of the final 
transition can be mitigated to some extent by assuming that the constant A^ in 
Eq. (2) is very small. In this case the static thermal recovery term grows more 
slowly according to the fourth power of the equilibrium stress and this produces a 
smoother transition in the stress-strain curve and hysteresis loops. 

The theory does not include 0 terms in the state variable evolution equations, 
so that the state variables cannot change during nonisothermal "elastic" excursions. 
If a theoretical yield surface is constructed from the model the surface can expand 
and translate in stress space but no rotation of the surface can occur. 

3.8 Cernocky and Krempl’s Theory 

This theory (Ref. 20) is written in an integral form in Appendix 10. In this 
particular form the theory is suitable for monotonic loading; for cyclic loading 
the "constants" E^, a, 3, K, 6 and ^ are updated discontinuously according to a 
set of rules. Some of these updating rules are presented in Ref. (21) and the 
others are currently under development. 
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Both an equilibrium stress and a drag stress K are utilized in the theory. 
The updating of E^., a and 6 in the equilibrium stress and the updating of the drag 
stress K with cyclic deformation allows for hardening in both state variables. 
Explicit relations are used for the equilibrium stress and the drag stress, so that 
only one differential (or integral) equation has to be solved (Eq. (1) in Appendix 
10) rather than three. At large strains the state variables saturate to values 
which are independent of strain rate. Static thermal recovery terms have not yet 
been specified in the theoretical formulation. 

^Presumably, since the "constants*' in the equilibrium stress, and in tfhe drag 
stress itself, depend on the current temperature, the model permits the state 
variables to change under nonisothermal "elastic" excursions. 

A theoretical construction of the yield surface from the model shows that the 
surface can both expand and translate in stress space, but cannot rotate. 


3.9 Hart’s Theory 

Hart’s theory (Ref. 22) is summarized in Appendix 11. Two state variables are 
employed, one being the equilibrium stress s the other a scalar variable, 

a*, called the "hardness" by Hart. The "hardness", a*, may be called a secondary 
state variable since it serves only to modify the equilibrium stress state variable, 

, at large strain values (or after the imposition of cyclic loading) in the same 
manner as the secondary state variables in Lee & Zaverl’s constitutive theory 
(Appendix 6) . 

The equilibrium stress grows linearly with inelastic strain in the initial 
loading phase and saturates when the static thermal recovery term containing the 
"hardness" state variable balances the linear work hardening term (the interpreta- 
tion of the "hardness" term as being a static thermal recovery term is the author’s). 
Since dynamic recovery terms are not included in the theoretical formulation, the 
saturated value of the equilibrium stress depends on the strain rate. The linear 
work hardening growth of the equilibrium stress together with the rapid growth of 
the static thermal recovery term at large strain values produces the tri-linear 
stress-strain response characteristic of Miller’s theory and Krieg, Swearengen and 
Rohde’s theory. A demonstration of this tri-linear behavior may be found in 
Fig. 6 of Delph’s review paper (Ref. 23). Since the drag stress is assumed to be 
constant the theory can model only kinematic hardening. Due to the changes in 
the "hardness" state variable the saturated equilibrium stress can change with 
cumulative deformation during cycling in the manner depicted in Fig. 5. Here the 
peaks of the stress-strain hysteresis loop increase in amplitude, the "yield stress" 
decreases in value, and the tangent modulus increases as cycling continues. This 
behavior is just noticeable in Fig. 6 of Delph’s paper (Ref. 23). It would have 
been even more apparent if Delph had increased the strain amplitude in his theo- 
retical predictions. The fact that the drag stress is constant implies that iso- 
tropic hardening is absent in the theory. Hence, at high strain rates, where the 
saturated stress values in the peaks of hysteresis loops have the approximate 
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values, o ^ ± Ke the hardening is not apparent. The relative amount of 

hardening will increase at lower strain rates where the stress has the approximate 
form a ± . 

The theory does not include 0 terms in the state variable evolution equatiofis, 
so that the state variables cannot change during nonisothermal "elastic” excursions. 
If a theoretical yield surface is constructed from the model the surface cannot 
expand or rotate in stress-space since only kinematic hardening is permitted. 


3.10 Valanis’ Endochronic Theory 

The endochronic theory (Ref. 24), summarized in Appendix 12, has been used 
mainly in a rate- independent form. Strain rate effects with creep-plasticity inter- 
action are still being developed by Wu and Chen (Ref. 25) and the rate-dependent 
aspect of the theory is still open. In the theory proposed by Wu and Chen (Ref. 25) 
different material functions are used in the formulation to describe the behavior 
of the material under constant strain rate, creep and relaxation conditions. A 
unification of the three material functions into one general function is required 
before the theory can be used under arbitrary loading conditions at high temperatures 


3.11 Laflen and Stouffer’s Theory 

This theory (Ref. 26) has no provision for cyclic hardening or static thermal 
recovery. It is therefore not suitable for describing the constitutive behavior 
of materials which exhibit significant cyclic hardening/softening. In its integral 
form the theory requires the evaluation of the entire integral of the strain history 
from the initiation of loading to the current time. For each new choice of the 
current time the entire integral must be recomputed (Ref. 28). This will exact 
severe economic and time penalties if it is used in a nonlinear finite element 
structural analysis code. 

Three theories were selected for further evaluation based on the constitutive 
review in this section. These are: (a) Walker’s functional theory; (b) Miller’s 

theory; and (c) Krieg, Swear engen and Rohde’s theory. However, the remaining 
theories are undergoing active development and it may be that in several years it 
will be possible to construct a hybrid theory which contains the best elements 
from each theory. 
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4.0 TASK II - CONSTITUTIVE TESTING 


4.1 Choice of Material 

Jet engine combustor liners are constructed from Hastelloy-X material. This 
material is a fine grained nickel base alloy which has relatively low strength but 
high ductility at elevated temperatures. The material constants required by the 
viscoplastic constitutive theories must be determined at various temperatures up 
to 2000°F (1093°C) in order to predict the structural response of combustor liners 
under engine operating conditions. 


4.2 Choice of Testing Machine 

The tests described in this section were conducted on uniaxial bar specimens 
using an INSTRON servo-hydraulic closed-loop testing machine capable of maintaining 
load or strain control with the inclusion of hold times. Specimens are maintained 
at a given temperature (to within 2°C (4°F)) in an oven heated by resistive elements 
which surround the specimen. Temperature control is maintained by a closed-loop 
system with thermocouples located at the ends of the cylindrical specimen. ASTM 
Class B-1 extensometry is used for the axial strain measurement. Load and total 
strain for this system can be controlled to within one percent. At + 0.6 percent 
strain amplitude, a strain rate of 3.87 x lO""^ sec ^ corresponds to traversing the 
hysteresis loop about ten times per minute. Since the temperature of the specimen 
is controlled by thermocouples placed at its ends, it is possible that the 
dissipation of inelastic work into heat at these strain rates could raise the 
temperature at the center of the specimen above that registered by the thermocouples 
located at its ends. A computation in Appendix 13 shows that at 982°C (1800°F) the 
rise in temperature of the center of the specimen above the temperature at its ends 
is only 0.8°C (1.4®F), so that the effects of inelastic dissipation of work into 
heat can be neglected in the data reduction of the Hastelloy-X tests. 


4.3 Uniaxial Hysteresis Tests 

Uniaxial bar specimens of Hastelloy-X (see Fig. 10) were tested at temperatures 
of 427®C (800°F), 537°C (lOOO^F) , 648°C (1200°F) , 760°C (1400^F) , 871^C (1600^F) and 
982^C (1800®F). 

This material is a fine grained nickel base alloy which has relatively low 
strength but high ductility at elevated temperatures. At each temperature the 
steady state hysteresis loops were determined under fully reversed strain controlled 
conditions. The input strain history consisted of a triangular wave function from 
a function generator with a strain amplitude of + 0*^ percent at 760°C (1400°F) , 
871^^0 (1600®F) , 982®C (1800°F) and a strain amplitude of + 1.0 percent at 427°C 
(800°F), 537°C (lOOO^F) and 648^C (1200°F) . 


25 



Hysteresis loops were obtained at strain rates of 1.25 x 10 1.11 x 10“ 


1.11 X 10“^, 

9 7 

3.70 x 10"^, 3.66 X 10“"^ and 3.87 x 10 sec“^ at each temperature. The material 
exhibits hysteresis loops which are asymmetric with respect to the strain axis at 
all temperatures, with the material being stronger in compression than in tension, 
and with the asymmetry increasing with temperature. However, the hysteresis loops 
are symmetric with respect to an axis which lies below the zero stress point. ^ 
Experimental results for Hastelloy-X are shown in Figs. 11-16. 


At 427®C (SOO^'F) and 537''C (1000°F) an inverse strain rate sensitivity 
(Figs. 11 and 12) was observed in which the peak stresses of the hysteresis loops 
decreased with increasing strain rate. Apparently, the change from normal to 
inverse strain rate sensitivity takes place between 537®C (1000°F) and 648^^0 
(1200°F) . It is observed that as the temperature decreases from 982°C (1800°F) 
to 648°C (1200°F) the "yield stress" of the material increases, so that the 
hysteresis loops at 648°C (1200''F) show a greater stress amplitude, for a given 
strain amplitude, than the hysteresis loops at 760®C (1400^F) , 871®C (1600®F) 
and 982®C (1800°F) . However, the "yield stress" decreases when the temperature 
decreases from 537"'C (1000®C) to 427®C (800''F) , and the stress amplitude at 537°C 
(1000°F) is larger than the stress amplitude at 427°C (800°F) for the Hh 1 percent 
strain amplitude hysteresis loops shown in Figs. 11 and 12. 


In their present forms the constitutive theories under consideration cannot 
predict the inverse strain rate sensitivity observed in Hastelloy-X at 427°C 
(800°F) an 537°C (1000®F) . However, since the strain rate sensitivity at these 
temperatures is relatively small, the material constants have been chosen to make 
the constitutive theories relatively rate insensitive at these temperatures. 


Miller (Ref. 26) has included inverse rate sensitivity in a recent version 
of his theory by including another state variable to account for the effects of 
solute hardening. Effectively, this changes the drag stress of the material, 
so tha^t at certain temperatures the drag stress decreases rapidly with increasing 
strain rates, and the material exhibits inverse rate sensitivity. In Section 3.1 
a constitutive formulation was proposed in which the response to infinitely fast 
loading conditions is given by a = E , where 5 is the instantaneous stress state 
variable whose growth law is prescribed. Inverse rate sensitivity may then be 
incorporated into the theory by assuming that the magnitude of 5 is smaller than 
the magnitude of the equilibrium stress 0 , under monotonic loading conditions at 
constant strain rate. Both Miller’s method and the method here proposed require 
the abandonment of an infinitely fast elastic response in order to predict inverse 
rate sensitivity. 
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4.4 Creep and Relaxation Tests 


Creep tests were performed at each temperature by starting at various points ? 
on a steady state hysteresis loop carried out under fully reversed strain controlled 
conditions at a strain rate of 3.66 x 10”^ sec”^ at + 0.6 percent strain amplitude. 

At a given point on the hysteresis loop the load was held at a constant value ahid 
the resulting creep behavior was recorded on a strip chart. After each creep » 

test the material was cycled at a strain rate of 3.66 x 10”"^ sec”^ until the 
hysteresis loop stabilized before proceeding to the next creep test. Asymmetry 
in the creep response was noted at each temperature, with the material being » 

stronger in compression than in tension, in conformity with the uniaxial hysteresis ' 
loop results. Figures 17 and 18 show the creep results at 871°C (1600°F) and 982°C 
(1800°F) , in which the strain at the beginning of each creep test has been normalized 
to zero to facilitate comparison with the theoretical predictions. 

In a steady state hysteresis loop executed at a constant strain rate, the 
equilibrium stress will saturate to a constant value at large strain amplitudes. 

This value can be determined from the preceding creep tests by locating the points 
on the unloading branches of the hysteresis loop at which the initial creep rate 
is zero. It should be emphasized that this point is difficult to define in an 
accurate manner due to the fact that the initial creep rate is small over a 
considerable portion of the unloading branches. A more accurate determination of 
this point can be made by observing the point at which the initial relaxation rate 
is zero. However, the function generator on the present experimental equipment 
would not allow a strain hold on the unloading portion of the hysteresis loop, 
so that the saturated equilibrium stress at a given strain rate and temperature was 
determined from the creep results. This saturated value of the equilibrium stress 
is used in the determination of the material constants. 

Combustor liners are constructed from Hastelloy-X sheet material. In order to 
deter;mine if the material constants developed from the Hastelloy-X bar specimens 
are appropriate for Hastelloy-X sheet material, creep tests were performed on 
Hastelloy-X sheet specimens at 871®C (1600°F) and 982°C (1800°F) . The sheet 
specimen is shoxm in Fig. 19. It was initially Intended to compare the monotonic 
uniaxial tensile response of bar and sheet specimens at these temperatures, but the 
small area of the sheet specimens precludes this comparison. During the initial 
heat-up of the sheet specimens from room temperature to 871®C (1600°F) and 982°C 
(1800^F) , it was necessary to keep the sheet specimens under a small tensile load 
to make sure that buckling did not occur. At these elevated temperatures the 
sheet specimens creep under the small tensile load before the final temperature 
is attained. Hence, a comparison of the uniaxial tensile stress-strain response 
was abandoned. Instead, creep tests were repeated at 871°C (1600°F) and 982°C 
(1800®F) to test the variation of secondary creep rates from specimen to specimen. 
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The sheet specimens were loaded from zero to full load in three seconds under 
load control on a Data-Trak system and allowed to creep at constant load at 871^C 
(ISOO'^F) and 982°C (1800°F) . Creep tests on the Hastelloy-X bar specimens were 
initiated at constant load from various points on a steady state hysteresis loop. 
Since the initial conditions in the bar and sheet specimens are different it is 
not appropriate to compare the primary creep behavior. However, under conditions 
of steady state secondary creep the initial conditions are not of concern, and an 
examination of steady state secondary creep rates will provide a valid comparison 
between bar and sheet creep behavior at 871®C (1600°F) and 982°C (1800®F) . 

The test results are shown in Tables 3 and 4. In interpreting these results 

it shbuld be kept in mind that due to the dependence of the creep rate on the 4th 

and 5th powers of the applied stress at 982®C (1800°F) and 871"^C (1600°F) , viz. 
c = Aa^ where n = 4 or 5, a small error in the stress (load) measurement can cause 
large errors in the creep rate. If the experiments were turned around and the 
saturated stress was measured due to a constant applied strain rate - in a 
monotonic tensile test for example - the agreement between theory and test would 
improve. Most of the secondary creep rates obtained with the sheet and bar 
specimens agree to within a factor of 1.5. Exceptions occur in the sheet specimen 

at a ~ 74 MPa (10700 psi) at 982°C (1800°F) and at a = 148 MPa (21500 psi) at 

871°C (1600°F) . The secondary creep rates in the sheet specimens were obtained 
from strip chart records. Variations in secondary creep rates by a factor of 1.5 can 
easily be obtained depending on which part of the chart the record is taken from. 

This figure provides an indication of the accuracy in reading from one record of 
one specimen. 

From the creep tests at 871®C (1600°F) and 982°C (1800°F) one can say that 
the material constants determined from Hastelloy-X bar specimens are suitable 
for representing the constitutive behavior of Hastelloy-X sheet material. This 
is based on the experimental agreement in strain rates at a given stress level 
obtained^ in creep and hysteresis testing. 

Relaxation behavior at 871°C (1600®F) and 982°C (1800°F) was determined by 
holding the strain at a zero value on the compressive portions of the steady state 
hysteresis loops carried out at 3.87 x 10”^, 3.66 x 10~'^, and 3.70 x 10"^ sec with 
a strain amplitude of + 0.6 percent. Representative results obtained at 982®C 
(1800°F) are shown in Figs. 20-22. 


4.5 Thermomechanical Test 

Work performed under Contract NAS3-21836 with the NASA-Lewis Research Center 
on a simulated combustor liner rig has shown that after about one thousand (1000) 
operation cycles (Ref. 8) fatigue cracks initiate at the louver lip (point A in 
Fig. 1) and propagate axially towards the seam weld (point C in Fig. 1). In order 
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to predict the number of cycles to crack initiation the life prediction algorithms 
require a knowledge of the stress-strain hysteresis behavior at the fatigue critical 
location at point A. Since strain gages will not operate cyclically at 982°C (1800°F)^ 
the experimental stress-strain hysteresis behavior at point A cannot be measured. 
However, the stress field is one dimensional at point A, and consists of a uniaxial 
component in the circumferential, or hoop, direction. Moreover, because the ^ 
combustor liner is a thermally driven structure, the total strain at the fatigue j 

critical location is zero and the mechanical strain history at this point can be 
determined from the temperature history 0(5) as a (5) [90(5) /95] 

Thermocouples on the louver lip provide the experimental temperature history. Thus, ^ 
both the mechanical strain history and the temperature history at the fatigue 
critical location, where the stress is uniaxial, are known. This temperature vs. 
mechanical strain history can then be applied cyclically to a uniaxial cylindrical 
tube specimen of Hastelloy-X to determine its stress vs. mechanical strain hysteresis 
behavior. Such a test is known as a "faithful cycle" test because it is meant to 
"faithfully" represent the "true" thermomechanical environment at a critical fatigue 
location in a combustor liner. 

A thermomechanical test was performed on Hastelloy-X with a standard closed 
loop servohydraulic test machine using low frequency (10 kHz) induction heating 
and compressed air for temperature control on the specimen (Ref. 27). Specimen 
temperature measurement provided by a radiation pyrometer, in conjunction with 
independently computer controlled preprogrammed mechanical strain and temperature 
histories, allow virtually any complex cycle to be applied to the specimen. The 
system is capable of cycling at positive or negative mean strain, with either 
stress or strain control and hold times within each cycle. The tubular specimen 
used in the thermomechanical test is shown in Fig. 23. Axial strain is measured 
from the machined internal ridges utilizing a linear variable differential 
transformer and a quartz internal extensometer . Load and total strain for this 
system can be controlled to within 1 percent and temperature to within 2°C (4°F) . 

Figure 24 shows the experimental behavior of Hastelloy-X in a thermomechanical 
cycle. The temperature history consists of a sine wave with a forty one (41) 
second hold period at the maximum temperature of 975°C (1750®F) as shown in 
Fig. 25, the lowest temperature in the cycle being 504°C (940°F) . The mechanical 
strain history is depicted in Fig. 26 and a corresponding cross-plot of the mechanical 
strain vs. temperature cycle is shown in Fig. 27. Axial mechanical strain in the 
specimen varied between -0.1 percent and -0.445 percent. 

It should be noted that in the thermomechanical test the actual input temperature 
history consists o^ a sine wave with a hold period at the maximum temperature point 
of the cycle. In the actual liner rig experiment the temperature history, measured 
by thermocouples located at the fatigue critical location (see Fig. 68), is not 
a sine wave with a hold period. To account for this slight difference in 
temperature history the strain rate in the thermomechanical test is adjusted so that 
the cross-plot of mechanical strain vs. temperature is the same in the thermo- 
mechanical test and in the computations produced by a structural analysis of the 
combustor liner rig. 
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4.6 Combustor Liner Rig Test 

A simulated combustor liner has been tested for fatigue life in a specially ^ 

'constructed experimental rig. This work was performed under Contract NAS3-21836 
with the NASA-Lewis Research Center and is reported in Ref. 8. Cyclic engine 
level temperatures and thermal gradients are generated on simulated combustor ^ 
liners by a 250 kW, 450 kHz induction heater used in conjunction with controlled ’ 

cooling air temperature and flow rate. A schematic drawing of the rig depicting 
the principal components is shown in Fig. 28. 

Cooling air is supplied to the lower plenum of the rig from a non-vitiated ' 

upstream air heater at temperatures up to 538°C (1000°F) • The air in the lower 
plenum is then directed through flow straightener plates to the upper flow divider 
plate which permits ratioing of the shroud-side airflow to the hot-side airflow. 
Shroud-side air provides the total cooling air supplied to the specimen which is 
discharged to the hot-side flow annulus, and then, together with the metered hot- 
side flow, ejected to ambient temperature. 

The specimen is inductively heated by power supplied from a 250 kW, 450 kHz 
induction heater. High frequency (450 kHz) was chosen to minimize depth of 
penetration of the induced current (heat) in the test specimen to best simulate 
the thermal surface loadings (radiative/conductive) experienced by combustor liners 
in engine operation. 

The induction coil is supported on a framework of glass-bonded mica which is 
transparent to the induction field. A 48.3 cm (19 in.) diameter quartz cylinder 
(also transparent to the field) is positioned between the coil and the specimen 
to form the inner boundary of the hot-side flow annulus. The quartz cylinder is 
retained in position by a 45.4 kg (100 lb) quartz cover. A cylindrical test 
specimen of Hastelloy-X, which simulates a combustor outer liner and incorporates 
five (5) sheet metal louvers, was selected for fatigue testing (Ref. 8). Prior 
to testing, the cylindrical louvered test specimen, transition pieces and top 
shroud, are assembled as a unit. This provides a permanent instrumentation 
installation and facilitates frequent inspection during testing. The assembled 
louvered test specimen and cover are shown in Fig. 29a. Fine wire chrome-alumel 
thermocouples are used to determine temperature distributions during testing. 

Typical thermocouple locations are shown in Fig. 29b. 

The test program consisted of a ninety (90) second thermal cycle in which the 
rig cooling air was maintained at a constant temperature and flow rate (of 504''C - 
2.5 kg/sec (940°F - 5.5 Ib/sec) while cycling power from the induction heater. 

A representative transient and steady state temperature distribution measured on 
the center (third) louver is shown in Fig. 26. The cycle consists of a twenty (20) 
second transient from an isothermal distribution of temperature at 504°C (940^'F) 
to a maximum temperature distribution, a forty (40) second steady state portion, and 
a cool down back to the original isothermal condition. 
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During the testing, a maximum circumferential temperature difference of 44 
(80'^F) was recorded at any time point in the cycle. For purposes of thermal 
and structural analyses, measurements were averaged to provide a circumferentially 
uniform temperature profile, 

A comparison of a theoretical transient heat transfer analysis with the .4 
transient temperature distribution determined by the thermocouples may be found 
in Ref, 8 . This thermal analysis was used to produce a thermal increment file 
for use in the structural analysis of the combustor liner rig specimen described 
in Section 6.7. Thermal histories measured by thermocouples on the specimen are 
shown in Fig. 30, Note that the thermal increment file used to ’’drive^* the three 
dimensional combustor liner analysis corresponds to the actual temperature in the 
combustor liner rig and is not a sine wave approximation to the temperature history. 
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TABLE 3 


Temperature 982°C (1800°F) 




Stress (psi) (MPa) 

Secondary Creep Rate 
(sec“^) for Sheet Material 

Secondary Creep Rate 
(sec*“^) for Bar Material 

28800 199 

Spec #1 2.5 X 10"^ 

- 

14300 99 

Spec #2 9.8 X lO"'^ 

9.0 X 10"'^ 

10700 74 

Spec #3 3.9 X 10“^ 



Spec #4 1.1 X 10"^ 

1.0 X lO”"'"'^ 

7200 50 

Spec #5 1.5 X 10"'^ 



Spec #6 8.3 X 10”^ 

9.0 X 10"^ 


TABLE 4 


Temperature 871°C (1600“F) 


Stress (psi) (MPa) , 

Secondary Creep Rate 
(sec~^) for Sheet Material 

Secondary Creep Rate 
(sec~^) for Bar Material 

21500 148 

Spec #7 6.1 X 10"^ 



Spec #8 1.1 X lO"'^ 

1.3 X 10“^ 

14200 98 

Spec # 9 1.6 X 10“^ 

c 


Spec #10 1.6 X 10"^ 

1.4 X 10“ 

7150 49 

Spec #11 2.2 X 10“^ 



Spec #12 2.4 X 10"^ 

1.5 X 10"^ 
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5.0 TASK III - INCORPORATION OF CONSTITUTIVE 
THEORIES INTO NONLINEAR FINITE ELEMENT CODE 


f 


5.1 Description of The MARC Program 

s 

The viscoplastic theories of Walker, Miller and Krieg, Swearengen and Rohde . 

have been incorporated into the MARC general purpose, nonlinear, finite element 
program, (cf. Ref. 1). This program, in common with other general purpose finite element 
programs, has been developed expressly for nonlinear structural analysis. These 
programs involve sophisticated computational algorithms and advanced finite ele- ^ 
meat formulations, yet rely on constitutive models whose applicability to the hot 
section component environment is questionable. 

The constitutive theories commonly in use in such finite element programs 
artificially partition the total strain into rate-independent plastic and rate- 
dependent creep components. These computer programs are then used to approximate 
the time-dependent inelastic behavior of a structural component by alternately 
applying increments of creep and plasticity. This artificial partitioning of the 
total strain into components which are governed by separate constitutive relation- 
ships has no underlying physical basis and does not take account of the fact that 
creep and plasticity should interact within the constitutive equation. 

A number of viscoplastic constitutive theories, in which creep and plasticity 
are united, have recently been proposed in the literature and were reviewed in 
Section 3.0. Scant experimental data exists to determine the material constants 
required by such theories. In addition, such theories comprise a system of "stiff” 
differential equations whose lack of stability has been a deterrent to their in- 
corporation into large general purpose finite element programs. 

In this contract the viscoplastic constitutive theories were incorporated into 
the MARC program by means of an initial stress technique. All of the material 
nonlinearity in the constitutive equations is incorporated into an initial load vec- 
tor and treated as a pseudo body force in the finite element equilibrium equations. 
Because the viscoplastic constitutive theories form a "stiff" system of differential 
equations, it is necessary to form the incremental constitutive equation appropriate 
to the finite element load increment by means of a sub increment technique. 

In the subincrement technique the finite element load increment is split into 
a number of equal subincrements and the viscoplastic constitutive theories are 
integrated over the small subincrements to form an accurate representation of the 
incremental constitutive equation over the finite element load increment. Integra- 
tion over each subincrement was accomplished by means of an explicit Euler forward 
difference method. Provided the subincrements are sufficiently small (so that the 
stability level of the forward difference method is not exceeded) , the technique 
has been found to work efficiently and accurately, even for large finite element 
load increments. However, the solution time required for a complete load increment 
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in a finite element program is linearly related to the number of sub increments used 
in formulating the initial load vector. There is, therefore, a considerable incen- 
tive to use as few sub increments as possible, consistent with the stability of the 
differential equations comprising the constitutive theory. 

A method for reducing the number of subincrements has been developed for the 
functional theory. The constitutive relations for this theory are written in an 
integral form and recursion relations have been developed to integrate the equations 
over a subincrement in an efficient, accurate and stable manner. Recursion relations 
for the constitutive theory are implicit and a Newton-Raphson technique is required 
within each subincrement to accurately integrate the constitutive theory over the 
sub increment . The integral relations are sufficiently accurate and stable to permit 
the use of only one sub increment per finite element load increment. In this case 
the functional theory is integrated over an increment equal in size to the finite 
element load increment. Inaccuracies may arise from large finite element load in- 
crements due to the fact that the strain increment vector and the temperature 
increment A0 are assumed to be constant over the finite element load increment. 

Hence the increments must be sufficiently small to enable the input history to be 
accurately modelled by piecewise constant values of AP^/At, As. /At and A©/ At, where 
AP^ is the applied incremental load vector. The fact that Ae^^At is assumed to be 
constant over the increment implies that proportional loading occurs over the incre- 
ment. Evidently the load increment must be small enough for this to be a valid 
approximation. 

The principle of virtual work may be used to generate the MARC nonlinear 
equilibrium equations governing the incremental response of the structure to an 
increment of load. In evaluating the nonlinear structural response of a component, 
the program assumes that the load history is divided into a number of incrementally 
applied load steps. Each load step is sequentially analyzed as a linear matrix 
problem using an appropriate stiffness matrix and load vector. Although each load 
step uses linear matrix methods to solve the incremental equilibrium equations, the 
incremental equilibrium equations themselves are nonlinear since the load vector 
will depend on the displacement increment obtained in the solution of the incre- 
mental equilibrium equations. 

The principle of virtual work may be written, for applied external point loads 
Pi, or displacements u^, in the form 

2^^ SeT cTj dv^Sujpj , (5.1) 

where the integral extends over the volume, V, of each finite element and the sum- 
mation sign extends to all the elements in the structure. 
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In Eq. (5.1) the virtual displacement vector is related to the virtual 
strain vector 6e^ through the relationship 


S€j=BjjSuj or S€pSuj bJ-, 


(5.2) 


where B.. is the strain-displacement matrix and the superscript T denotes trans- 
position‘s Since is an arbitrary virtual displacement vector, Eqs. (5,1) and 
(5.2) may be written in the form 

(5.3) 

This relation expresses the equilibrium of the structure when the applied load 
vector is and the stress vector is a^. If an incremental load AP^ is applied 
to the structure and the stress vector changes to + Aa^, the relation expressing 
the equilibrium of the structure at the end of the incremental load application 
may be written as 


Z j^B|j(o-j+A<rj)dV = P, +APj • (5.4) 

Hence, the relation expressing the equilibrium of the structure during the applica- 
tion of the incremental load vector AP^ is obtained from Eqs. (5.3) and (5.4) by 
subtraction in the form 


Z/bJ; A(t, dV = AP; • (5.5) 

The MARC code allows the user to implement very general constitutive relation- 
ships into the program by means of the user subroutine HYPELA. Within this sub- 
routine the user must specify the values of the elasticity matrix and the 
inelastic stress increment vector ACj^ in the incremental vector constitutive rela- 
tionship 


Ao-j = Dj|(A€. -8jO A®) - A^. • 


(5.6) 


The inelastic stress 
viscoplastic constitutive 


increment vector A?, is computed in HYPELA from the 
relationships summarized in the Appendices. 
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In Eq. (5.6) a denotes the coefficient of thermal expansion and 6. is the 
vector Kronecker delta symbol. 


(5.7) 

For the class of nonlinear visoplastic constitutive relationships under con- 
sideration in this contract, the incremental inelastic stress vector depends 
in a highly nonlinear manner on the incremental strain vector Ae^^. Since Ac^ = 

Bj[j Auj , the incremental inelastic stress vector Ac^ depends in a highly nonlinear 
manner on the nodal displacement vector increment j > so that A^^ = A^^ (Auj). 

Substitution of Eq.. (5.6) into (5.5) produces the incremental equilibrium 
equations for MARC in the form, 

I Kjj Auj = APj + ARj + I^B^'j ACjdV + ZjfBj’jSjaA© dV , (5.8) 

where is the elemental elastic stiffness matrix defined by the relation 

= Dkf dV. (5.9) 

V ^ 

The vector AR^ is the residual load correction vector or out-of-equilibrium 
force vector from the preceding load increment, 

ARpPj-Zj* BjjCTjdV , (5.10) 

which is added to the current increment in order to restore the structure to 
equilibrium. The nonlinearity in the incremental equilibrium relationship, de- 
fined in Eq. (5.8), arises because the inelastic stress increment vector A^^ 
depends nonlinearly on the displacement increment vector Auj . Values of and 
ACj appropriate to the current incremental load step are returned to the main pro- 
gram by subroutine HYPELA and the incremental equilibrium relations in Eq. (5.8) 
are solved by successive iterations. 

The solution of the incremental equilibrium equations in (5.8) is accomplished 
within the MARC code by the following algorithm. At the start of the increment 
the user subroutine HYPELA is entered to determine the elasticity matrix D^j and 
the incremental inelastic stress vector A^ ^ . On entry to the subroutine the input 
consists of the strain increment vector Ae^, the temperature increment A0, the 
time increment At over which the incremental external load vector AP^ is applied 




'I if 0< j< 3 
"Oif3<i<6 . 
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to the structure, and the values of the stress, strain, temperature and viscoplastic 
state variables at the beginning of the increment. Since the incremental strain 
vector, Ae. = B..Au., can only be accurately determined after the solution to the ^ 
incremental equilibrium relationship in Eq. (5.8) has yielded the correct incremental 
solution vector Au. , the strain increment vector Ae . initially used to generate the 
inelastic stress increment vector A^. must be estimated* The initial estimate for 
is assumed to be the value obtained for Ae^ in the preceding increment* On exit’ 
from subroutine HYPELA the elasticity matrix and the estimated inelastic stress 
increment vector Ac. are passed into the main program* After the values of and 
Ac* are obtained for each integration point in the structure, the incremental equi- ^ 
librium relationship in Eq* (5*8) is assembled and solved for the incremental node 
displacement vector Au.* The incremental strain vector, Ae^^ ” ^ij then com- 

puted and compared wit^ the initial guess for Ae^ used to generate the inelastic 
incremental stress vector ACj . If this incremental strain vector is equal, within 
a user specified tolerance, to the incremental strain vector used to compute ACj 
in the assembly phase, the solution is assumed to have converged* Otherwise the 
updated strain increment vector, obtained from the solution of the equilibrium re- 
lations in Eq. (5.8), is passed into subroutine HYPELA, a new vector, ACj> is com- 
puted, and the equilibrium equations resolved to yield an improved value of Au^ and 
Ae . . The process is repeated until the value of the vector Ae^ on the assembly 
phase is equal, within a user specified tolerance, to the value of the vector Ae^^ 
on the solution phase* After convergence is achieved, the temperature, stress 
vector, strain vector and viscoplastic state variables are updated by adding the 
incremental values generated during the current increment to the values of these 
variables at the beginning of the increment. The program then passes on to the 
next load increment where the process is repeated* A flow chart of the iterative 
procedure required to implement the viscoplastic constitutive theories into the MARC 
program is shown in Fig. 31* 


5.2 Implementation of Functional Theory in MARC 


The recursive integration algorithm for the functional theory described in 
Appendices 14 and 15 has been introduced into the MARC code by means of the user 
subroutine HYPELA. This routine is called at each integration point in each 
element and supplies the elasticity matrix and the inelastic stress increment 


vector 




to the main program. 


The required header cards are: 

SUBROUTINE HYPELA ( D , G , E , DE , S , TEMP , DTEMP , NGENS , N , NN , KC , MAT , NDI , 
INSHEAR) 

DIMENSION D (NGENS, NGENS) , G (NGENS) ,E (NGENS) ,DE(NGENS) ,S(NGENS) 
DIMENSION TEMP (1) , DTEMP (1) 


FORTRAN CODING IN APPENDIX 17 


RETURN 


END 


37 



where 


through 


through 


D(NGENS,NGENS) 

is the elasticity matrix D^. defined in this 
subroutine (output argument; , 

G(NGENS) 

is the inelastic stress increment vector A?. ^ 

defined in this subroutine (output argument^, 

E(NGENS) 

is the mechanical strain ~ j ^a(0 [30(0 /9^ 

at the beginning of the increment (input argument) , 

DE(NGENS) 

is the increment of mechanical strain - 6^aA0 

(input argument), 

S(NGENS) 

is the stress at the beginning of the increment 

(input argument) , 

TEMP(l) 

is the temperature 0 at the beginning of the 
increment (input argument) , 

TEMP(2) 

is the time t at the beginning of the increment 
(input argument) , 

TEMP (3) 

is the cumulative inelastic strain R at the 
beginning of the increment (input argument) , 

TEMP (4) 
TEMP (9) 

are the values of the inelastic strains, Cj^ 
through Cg, at the beginning of the increment 
(input argument) , 

TEMP (10) 
TEMP (15) 

are the values of the equilibrium stresses 
through at the beginning of the increment 

(input argument) , 

TEMP (16) 

is the drag stress K at the beginning of the 
increment (input argument) . On exit from the 
subroutine TEMP (16) contains the drag stress K 
at the end of the increment (output argument) , 

DTEMP(l) 

is the temperature increment A0 (input argument) , 

DTEMP(2) 

is the time increment At (input argument) , 

DTEMP(3) 

is the increment of cumulative inelastic strain 
AR (output argument) , 
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DTEMP(4) 
through DTEMP(9) 

DTEMP(IO) 
through DTEMP(15) 

DTEMP(16) 

NGENS 


N 

NN 

KC 

MAT 

NDI 


NSHEAR 


are the incremental values of the inelastic 
strain Acj^ through Ac^ (output argument) , 


are the incremental values of the equilibrium 
stress Afi^ through Afi^ (output argument) , 


is output as zero, since the drag stress K is 
updated in TEMP(16), 


is the size of the D^. matrix (NGENS = 3 for plane 
stress problems, NGENS = 4 for plane strain and 
axisymmetric problems, NGENS = 6 for three dimen- 
sional problems) (input argument) , 

is the finite element number (input argument) , 

is the integration point number (input argument) , 

is the layer number in shell or beam problems 
(input argument) , 

is the material identifier (input argument) , 

is the number of direct stress components (NDI = 2 
for plane stress problems, NDI = 3 for plane strain, 
axisymmetric and three dimensional problems) (input 
argument) , 

is the number of shear components (NSHEAR = 1 for 
plane stress, plane strain and axisymmetric problems, 
NSHEAR = 3 for three dimensional problems) (input 
argument) . 




5.3 Notes on Subroutine HYPELA 

The values of D^^ and A?^ in the incremental constitutive relation. 


Ao-.. = Di|(A€. -SjO A®) - A^. . 


(5.11) 


are obtained by a sub increment method. Incremental values of the variables 
A0, At and (Ae. - 6^aA0) for the current finite element load increment are split 
into N equal values, and the constitutive equations are integrated over the N sub- 
increments to provide accurate values of D^j and hi;^. Because the recursive scheme 
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is accurate and stable, even for large load increments, it is possible to use only 
one sub increment per MARC increment. Figure 32 shows the MARC finite element pre- 
dictions for an axis 3 rmmetric problem is which a bar of Hastelloy-X material is 
strained at a uniform strain rate in tension. The rounded curve was obtained by 
loading to 0.64 percent strain in 64 load increments with one subincrement per MARC 
load increment. This corresponds to integrating the recursive algorithm directly 
over the same step size as the finite element step size. Each of the curves marked 
8, 4, 2 and 1, respectively, was obtained by using only one sub increment per MARC 
load increment and loading to 0.64 percent strain in 8, 4, 2 and 1 finite element 
load increments, respectively. It is evident that the recursive scheme is stable 
and accurate even when loading to 0.64 percent strain directly in only one MARC 
increment. 

The accuracy and stability of the integration algorithm will enable large in- 
crements to be used in a nonlinear finite element analysis with a concomitant re- 
duction in computer run time. Although the recursive integration algorithm is 
implicit, requiring Newton-Raphson iterations within each subincrement, the number 
of iterations required is not large. In the numerical algorithm the initial guess 
for the inelastic strain increment in the first subincrement is taken as the devia- 
•toric mechanical strain increment. If the material is in the elastic region this 
is not a good initial guess. However, if the response is elastic, then after two 
iterations it is found that the magnitude of the inelastic strain increment will 
generally be less than one percent of the mechanical strain increment and the itera- 
tions are stopped. If the material is in the inelastic region the initial guess 
for the inelastic strain increment is very well approximated by the deviatoric me- 
chanical strain increment and one iteration is sufficient to produce an accurate 
result. In successive suId increments the initial guess for the inelastic strain 
increment is taken from the preceding sub increment . Although the computer time 
required to integrate the constitutive relations over one subincrement is longer 
than the time required to integrate the same equations by means of an Euler forward 
difference method (by about a factor of three or four in general) , the recursive 
method i? stable and accurate even when the constitutive equations are integrated 
over a large finite element increment with only one subincrement. Moreover, the 
forward difference method has a stability limit, which precludes the use of large 
sub increments. Backward difference and implicit integration methods could also be 
used to integrate the constitutive equations over a sub increment . Although these 
methods are stable they are not as accurate as the recursive integration algorithm. 

MARC solves the incremental equilibrium Eq. (5.8) by successive iteration. 

To see how the equilibrium equation iterations are converging one can pick the 
integration point, NPRIN, in element number, NELPR, at which the maximum amount 
of nonlinearity is expected. As subroutine HYPELA is entered on the assembly 
phase the routine prints out the strain increment vector and the stress incre- 
ment vector at integration point NPRIN in element number NELPR. After the 
equilibrium equations have been solved for the incremental displacement vector Auj^, 
subroutine HYPELA is again entered with Ae^^ = B..Au. and the incremental vectors 
are printed out on the assembly and solution phase if every successive iteration of 
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the equilibrium equations. In this way the convergence of the solution to the 
incremental equilibrium equations can be followed. If no printout of the incre- 
mental vectors is required, the variable IPR is set equal to zero in subroutine 
HYPELA. 

For each subincrement the integral recursive algorithm requires Newton-Raphson 
iterations. A test is made within the subroutine to determine if the solution 
is in the elastic region. If it is in the elastic region and two Newton-Eaphson 
iterations have been performed, the iterations are stopped and the subroutine goes 
on to the next sub increment . Another test is also made to see if the minimum number 
of iterations is met, or until the tolerance on the magnitude of the inelastic strain 
increment vector, 

viz. I "" ^^-ll^^^n-l| - where n is the iteration number, 

is met. If the foregoing tolerance on the inelastic strain increment vector is 
not met but the maximum number of Newton-Raphson iterations, MAXIT, is exceeded, 
the subroutine exits from the iteration loop and proceeds to the next sub increment . 

In order to use subroutine HYPELA, nine (9) constants must be defined in the 
subroutine, starting at card number forty-four (44). These constants are: 


NTP 

MAXIT 

MINIT 

NELPR 

IPR 


number of tabulated temperature points in the 
DATA statements, 

maximum number of Newton-Raphson iterations 
allowed , 

minimum number of Newton-Raphson iterations 
allowed , 

element number for printout of incremental stress 
and strain vectors, 



1 if stress-strain increment output is required, 

0 if stress-strain increment output is not required. 


NPRIN 


integration point number for printout of 
incremental stress and strain vectors. 


NSPLIT = number of subincrements per MARC increment, 

SFTEMP = stress free temperature; for isothermal calculations 

this is the isothermal temperature. 


TDIF 


temperature difference in the DATA statement in 
which the material constants are tabulated. 
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The material constants for the functional theory, viz. 

O 

A, li, a, n, m, n^, n^, n^, , K^, in Appendix 1, 

depend on temperature. For Hastelloy-X these constants have been determined at 
six temperatures ranging from 427 (800®F) to 982°C (ISOO'^F) in lll^C (200°F) ^ 
intervals. These material constants are entered into the MARC subroutine HYPELA 
in DATA statements. At a specified temperature, the material constants are found 
by a linear interpolation of the constants in the DATA statements. If the tempera- 
ture lies outside of the tabulated values, a linear extrapolation, based on the 
last two entries in the DATA statement table, is performed. The material constant \ 
n is not entered directly into the DATA statements. Instead, values of its 
reciprocal, n are tabulated. At a given temperature n~^ is then found by linear 
interpolation in the table and n is determined as its reciprocal. 

If the material constant n is determined by linear interpolation of its 
tabulated values, unreasonable results are obtained for the material response. 

Figure 33 shows the steady state hysteresis loops for Hastelloy-X at 704°C (1300®F) 
based on a linear interpolation of the tabulated material constant n at 648®C 
(1200®F) and 760^C (1400®F) , This figure shows that the prediction of the 704°C 
(1300°F) response, based on the linear interpolation of n, lies outside of the range 
of the 640°C (1200®F) and 760°C (1400°F) responses. The constants at 648"^C (1200°F) 
and 760^C (1400°F) were obtained from the experimental data and the theoretical 
predictions at these temperatures in Fig. 33 provide a reasonable correlation with 
the experimental data. 

The tabulated material constants are obtained from the asymptotic response of 
the material at strain values large enough for the stress to obtain a saturated 
limit. In order to obtain a realistic response at a given temperature, the material 
constants should be obtained by a linear interpolation of the asymptotic response 
of the material. If and denote the theoretical asymptotic stress values cor- 
responding to a strain rate e at temperatures and © 2 , where the strain value e is 
large enough for the equilibrium stress to reach its maximum rate-independent value, 

0. , then 

max 


cT|=nn,o^{®,)+K(®,)e 


l/n(®|) 


(5.12) 


and 

From these equations it is clear that in order to obtain a realistic asymptotic 

response for a at a temperature 0 between 0]_ and © 2 > n linear interpolation of 

9. , K and n”^ is appropriate. Numerical calculations show that a linear inter- 

max 

polation of n is inappropriate, and linear interpolation of n is therefore per- 
formed in the HYPELA subroutine. 
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5.4 Input Data Required by MARC 


A listing of the data cards required for the thermomechanical loading history 
described in Section 4.3 is given in Appendix 16. This problem consists of the 
thermomechanical loading of a bar of Hastelloy-X material in which the bar is 
subjected to the mechanical strain and temperature histories depicted in Figs. 25 - 
27. Since only mechanical strain is being imposed on the bar, this condition can 
be simulated by setting the coefficient of thermal expansion a equal to zero and 
imposing a given total strain history at the ends of the bar. A four element 
axisymmetric problem is depicted in Fig. 34. Nodes one, two and three are given 
incremental displacements which follow the total strain history of Fig. 26, whilst 
ea,qh integration point of each element is given incremental temperature variations 
which follow the temperature history of Figure 25. 

Two cards are required before the END card (Card 10 of Appendix 16) . These 

are: 


HYPOELAS 

STATE VARS 16 

A uniform temperature increment over the structure, together with an appro- 
priate time increment, can be specified with the following cards; 

THERMAL LOADS 

1 , 

5.0, 2.0 

BLANK CARD 

In the above sequence of cards, the first state variable increment of 5.0 refers 
to the uniform temperature increment A0 = 5°F over the structure. The second state 
variable increment of 2.0 refers to a time increment of At = 2 seconds. The re- 
mainder of the card and the following BLANK card set the remaining fourteen (14) 
state variable increments to zero. Since the STATE VARS card (Card number 7 of 
Appendix 16) defines sixteen (16) state variables, MARC expects this number as 
input. However, only the first two state variables, A0 and At, are required as 
input by HYPELA, so that the remaining (14) state variable increments are set to 
zero. The remaining (14) state variable increments are defined within subroutine 
HYPELA. 

If a nonuniform temperature over the structure is required (the usual case) , 
the temperature increments and time increment must be set in user subroutine CREDE. 
This can be accomplished with the following header cards: 

SUBROUTINE CREDE(DTDL,M,NSTRES ,NEQST,NSTATS) 

DIMENSION DTDL(NSTATS,NEQST,NSTRES) 

COMMON/ FAR/DUM,L 

N = (where N = number of integration points in element number M) 
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DO 2 II = 1, N 


DTDL(1,1,II) = temperature increment at integration point II 

DTDL(2,1,II) = time increment (can be made a function of load 

increment number L) 

2 CONTINUE 
RETURN 
END 
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6.0 TASK IV - COMPARISON AND EVALUATION OF CONSTITUTIVE THEORIES 


6.1 Determination of Material Constants 

Material constants for each of the three candidate theories were determined 
from the experimental hysteresis data presented in Figs. 11-16 and from the satu- 
rated value of the equilibrium stress at a constant strain rate magnitude of 
3.66 X 10“*^ sec"^ at a temperature of 982®C (1800°F) . Values of the constants 
for the functional theory, Miller’s theory, and Krieg, Swearengen and Rohde’s 
theory are presented in Tables 5-7. 

The constant A^ in Krieg, Swearengen and Rohde’s theory has been assumed to 
have the value A^ = 10““^^ for all temperatures. The thermal recovery term may 
then be approximated in the form 


recovery = 


A2 A3 





( 6 . 1 ) 


so that in uniaxial loading the recovery term is proportional to the fourth power 
of the equilibrium stress. If A 3 were larger in value, the exponential term would 
grow very rapidly once the equilibrium stress reached a certain value, as shown in 
Fig. 9, and the resulting stress-strain curves would exhibit a tri-linear behavior. 
The small value of A 3 and the resulting dependence of the static thermal recovery 
term on mitigates the severity of this tri-linear behavior. 

The constants (n^ + n ) in the functional theory, in Miller’s theory and 
A^^ in Krieg, Swearengen and Rohde’s theory determine the rapidity of the approach 
of the equilibrium stress to its saturated value. For large values of these con- 
stants the equilibrium stress grows so rapidly that it virtually saturates in the 
’’elastic" region and the stress-strain curve exhibits the ’’square’’ behavior de- 
picted in Fig. 35. For intermediate values the equilibrium stress approaches its 
saturated value less rapidly and the stress-strain curve exhibits a rounded appear- 
ance. For very small values of these constants, the equilibrium stress grows very 
slowly. The asymptotic stress value, when ~ 0, is given by the functional theory 
and by Krieg, Swearengen and Rohde’s theory as 

( 6 . 2 ) 


and by Miller’s theory as 




i/n ^2/3 


(6.3) 
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TABLE 5 


Material Constants for Functional Theory 


Material 

982'’C 

871°C 

760°C 

648°C 

537°C 

427°C 

Constant 

1800°F 

1600°F 

1400°F 

1200°F 

1000°F 

800°F 

A 

11.5E6 

15.4E6 

17.8E6 

18.1E6 

17.2E6 

17.8E6 

y 

4.9E6 

6.9E6 

8.4E6 

9.0E6 

9.0E6 

9.8E6 

^1 

59292 

91505 

251886 

95631 

75631 

50931 

^2 

0 

0 

0 

0 

0 

0 

n-1 

.233 

.195 

.244 

.079 

.059 

.059 

m 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 


0 

0 

0 

0 

0 

0 

^2 

1.0E6 

5.0E6 

2.0E7 

1.5E7 

6.0E7 

30.0E7 

ns 

312 

673 

1179 

781 

1000 

8000 


0 

0 

0 

0 

0 

0 

^5 

0 

0 

0 

0 

0 

0 

^6 

2.73E-3 

8.98E-4 

0 

0 

0 

0 

^7 

0 

0 

0 

0 

0 

0 

0 

-1200 

-1434 

-2000 

-2000 

0 

0 
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TABLE 6 


Material Constants for Miller's Theory 


K 

= 

8000 

0 

n 


1.598 

B 

= 

1.0293E14 

H - 

= 

1.0E7 

1 

A 


9.305E-4 

1 

»2 


100 

^2 

= 

50000 

A 

ss 

5.9425E-12 


Q* = 104600 

T = 1588°K 

m 

k = 1.9859 

0' = exp I - Q*/kT I for T > .6T^ 

( - 0 * / . 6T„. \ ) 

e- = exp j In (l + _ for T < .6T^ 

where T is the temperature in degrees Kelvin, The Lame constants X and y 
are as given in Table 5. 
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TABLE 7 


Material Constants for Krieg, Swearengen and Rohde’s Theory 


Material 

Constant 

982°C 

1800°F 

871’’C 
1600 '"F 

760°C 

1400°F 

648°C 

1200°F 

537‘’C 

1000°F 

427°C 

800°F 

A 

11.5E6 

15.4E6 

17.8E6 

18.1E6 

17.2E6 

17.8E6 

y 

4.9E6 

6.9E6 

8.4E6 

9.0E6 

9.0E6 

9.8E6 

K 

'‘O 

59292 

91505 

251886 

95631 

75631 

50931 

n““^ 

.223 

.195 

.244 

.079 

.059 

.059 


1.0E6 

5.0E6 

2.0E7 

1.5E7 

6.0E7 

30.0E7 

A 

2 

243 

14.96 

1.54 

. 66 

1.79E-3 

.59 

A 3 

l.OE-12 

l.OE-12 

l.OE-12 

l.OE-12 

l.OE-12 

l.OE-12 


0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 
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Thus, for very small values of 112 the theories predict a viscoelastic material 
behavior and the stress-strain curve again exhibits a relatively ^square” behavior. 
Numerical results for the functional theory are exhibited in Figs. 36 and 37 when 
n^^ = 0 and n^ = 0.8E7 and n 2 = 5.0E7, respectively. The ^square" behavior for 
large values on n^ is apparent in Fig. 37. In changing the value of U 2 from 0.8E7 
to 5.0E7, it is necessary to adjust the material constant n^ in the dynamic recovery 
term in order to obtain the same saturated value of the equilibrium stress. 


6.2 Strain Rate Behavior 

"Figures 38 - 58 depict the steady state hysteresis loops obtained with the func- 
tional theory, Miller’s theory, and Krieg, Swear engen and Rohde’s theory at tempera- 
tures ranging from 648°C (1200'^F) to 982°C (1800°F) in 38"^C (100°F) intervals. 

These loops represent the predicted behavior of the material for triangular strain 

o —'4 

histories in which the magnitudes of the strain rates are 3.9 x 10"" sec , 3.7 x 10 
sec"^, 3.7 X 10“^ sec”^, 1.1 x 10”^ sec”*^ and 1.2 x 10““^ sec“^ and in which the 
strain amplitude is + 0*6 percent. Tri-linear material behavior is observed in 
Miller’s theory in the temperature range extending from 648°C (1200°F) to 816"^C 
(1500°F) , whilst the functional theory and Krieg, Swearengen and Rohde’s theory pro- 
duce rounder stress-strain hysteresis loops. 

Figures 59 and 60 facilitate a comparison between the functional theory and 
the experimental hysteresis loops at 871°C (IbOO^^F) and 982°C (1800®F). The peak 
stresses in the loops are well predicted, but at the lower strain rates the theoret- 
ical predictions of the stress-strain behavior are too "square". This "squareness" 
at the lower strain rates is common to each of the theories and arises from the fact 
that the equilibrium stress saturates to a smaller value at low strain rates due 
to static thermal recovery. For a given hardening rate of the equilibrium stress, 
the saturated value is achieved more rapidly at low strain rates, and this produces 
the "square" stress-strain loops. The theoretical loops at low strain rates could 
be made rounder by assuming that the constants which govern the hardening of the 
equilibrium stress, viz. (n^^ + ^ 2 ^^ ^ 1 * decrease in magnitude in proportion 

to the decrease in the saturated value of the equilibrium stress. This is equiva- 
lent to the assumption that these hardening constants depend on the magnitude of 
the inelastic strain rate, 

Figures 61 and 62 show the results of a MARC axisymmetric finite element 
analysis of a cylindrical specimen of Hastelloy-X using the functional theory. 

This specimen was cycled under fully reversed strain controlled conditions with a 
strain amplitude of +^.6 percent at 982°C (1800°F) at a constant strain rate magni- 
tude of 1.1 X 10“-^ sec”^. The MARC plots represent loading the specimen to +0.6 per- 
cent strain, unloading and compressing the specimen to -0.6 percent strain, and 
finally reloading the specimen to +0.6 percent strain. In Fig. 61 sixty (60) MARC 
increments were used to load from 0 to +0.6 percent strain, one hundred and twenty 
(120) MARC increments to load from +0.6 to -0.6 percent strain, and one hundred and 
twenty (120) MARC increments to reload from -0.6 percent to +0.6 percent strain. 
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with one sub increment per MARC increment in subroutine HYPELA. The reason for 
using a large number of increments to traverse the hysteresis cycle was simply 
to obtain good definition of the stress-strain hysteresis loops in the plotted 
output from MARC. Figure 62 shows the same hysteresis loop in which three (3) 
MARC increments were used to load the specimen from 0 to +0.6 percent strain, 
six (6) MARC increments to unload from +0.6 percent to -0.6 percent strain, and 
six (6) MARC increments to reload from -0.6 percent to +0.6 percent strain. 

Thirty (30) subincrements were used per MARC increment in subroutine HYPELA, A 
comparison of Figs. 61 and 62 shows that even with the large MARC increments used 
in Fig, 62 the stress at the end of the MARC increment lies on the stress-strain 
hysteresis loop. Experience with the algorithm shows that even if only one sub- 
increment per MARC increment is used, the stress at the end of the MARC increment 
li'Cs on the stress-strain hysteresis loop. 


6.3 Creep and Relaxation Behavior 

A comparison of experimental and theoretical creep predictions using the 
functional theory is shown in Figs. 17 and 18. The experimental results were gen- 
erated by traversing a steady state hysteresis loop carried out under fully reversed 
strain controlled conditions at a strain amplitude of +0.6 percent at a constant 
strain rate magnitude of 3.7 x 10~^ sec”^ at temperatures of 871''C (1600 °F) and 
982°C (1800°F). At various points on both the tension going and compression going 
loading branches of the hysteresis loops, the stress was held at a constant value 
and the creep results were recorded on a strip chart. After each creep test was 
completed the material was cycled again at +0.6 percent strain amplitude at a 
strain rate magnitude t>f 3.7 x 10”^ sec“^. After a few cycles the material will 
cycle around the original steady state hysteresis loop and the effects of the pre- 
vious creep tests are wiped from the materials memory. 

The functional theory was integrated around the hysteresis loops at 871°C 
(1600°F) and 982°C (1800°F). At various points on the hysteresis loops the stress 
was held at a constant value and the resulting creep strain produced by the vis- 
coplastic theory was computed. The dots in Figs. 17 and 18 represent the predic- 
tions of the functional theory. Since the theoretical hysteresis loops are too 
"square” in comparison with the experimental curves, the initial strain at which 
creep commences differs in the experimental and theoretical creep curves. To avoid 
this discrepancy the theoretical and experimental creep curves in Figs. 17 and 18 
have been normalized to a zero initial strain. 

According to the theoretical formulation, steady state creep is indistinguish- 
able from the stress-strain behavior at large strain values on a hysteresis loop 
executed at a constant strain rate, since at large strain values the stress satu- 
rates to a constant value. This conclusion is supported by the close agreement 
between the theoretical and experimental creep behavior, since the material constants 
were determined from the hysteresis loops and not from the creep curves. A differ- 
ence in the compressive and tensile hysteresis behavior found in the stress-strain 
loops is also evident in the creep behavior, the material being stronger in 
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compression than in tension. The introduction of the constant, ^2, into the 
functional theory to account for the different tensile and compressive hysteresis 
behavior also accounts for the difference between the tensile and compressive 
creep behavior shown in Fig. 18. 

Figure 63 depicts the creep response obtained with Miller’s theory at 1600 °F. 
Since Miller’s theory predicts equal response in tension and compression, only ^ 
the results for the compressive calculations are given. At low stress levels the 
predictions are not as good as the functional theory predictions. However, a 
reasonable creep prediction is obtained at the highest compressive stress level 
of -197 MPa (-28,600 psi) . These creep predictions have been made at a temperature 
where Miller’s theory should be at its most accurate. The accuracy obtained at 
other temperatures will decrease because the material constants provide a best fit 
at 1600°F and give considerable scatter in attempting to correlate the behavior at 
other temperatures. 

The creep behavior obtained with Krieg, Swear engen and Rohde’s theory at 871° C 
(1600°F) and 982°C (1800°F) is very similar to the creep behavior of the functional 
theory in the case of tensile creep. It predicts equal tensile and compressive 
creep behavior. 

Stress relaxation data was obtained at temperatures varying from 648°C (1200°F) 

to 982°C (1800°F). In the tests stress relaxation was initiated at a zero strain 

value on the compressive loading branches of three steady state hysteresis loops 

— —1 —4 —1 

carried out at constant strain rate magnitudes of 3.9 x 10 sec , 3.7 x 10 sec 
and 3.7 x 10“*^ sec~^ with strain amplitudes of +0.6 percent. After each stress 
relaxation test the material was cycled several times around the steady state 
hysteresis loops to erase the preceding hysteresis and relaxation tests from the 
memory of the material. 

Since the relaxation data was not used in the determination of the material 
constants, very little of the experimental relaxation data has been compared with 
the theoretical predictions to date. 


In Figs. 20 - 22 the dots represent calculations of the stress relaxation 
obtained by integrating the functional theory around steady state hysteresis loops 
executed at strain rates of 3.9 x 10“^ sec“^, 3.7 x 10 ^ sec”^ and 3.7 x 10 ^ sec ^ 
under strain controlled conditions with strain amplitudes of +®.6 percent at 982°C 
(1800°F). At zero strain values on the compressive portions of these loops the 
strain was held at zero and the resulting stress relaxation produced by the theory 
was computed. The agreement between the theoretical predictions and the experimental 
relaxation tests show satisfactory agreement. 


6.4 Thermomechanical Behavior 

The thermomechanical loops predicted by the functional theory. Miller’s theory 
and Krieg, Swearengen and Rohde’s theory are displayed in Figs. 64-66. These 
hysteresis loops are produced by integrating the theories over the temperature and 
mechanical strain histories depicted in Figs. 25-27. 
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From Fig. 64 it is evident that the functional theory prediction exhibits a 
relatively small amount of stress relaxation at the maximum temperature of 954°C 
(1750°F) in Fig. 25. Note from Fig. 26 that the maximum temperature hold period 
does not correspond to the maximum compressive mechanical strain. In fact the 
strain is held constant on the unloading branch of the thermomechanical hysteresis 
loop. The hysteresis loops show a progressive overall stress relaxation in the 
positive stress direction and the loops stabilize after about three cycles in ^ 
conformity with experimental observations. 

Figure 65 shows the thermomechanical loop predicted by Miller’s theory. The 
theory predicts a large amount of stress relaxation at the maximum temperature hold 
period and the loop stabilizes in one cycle. In the tensile portion of the thermo- 
mechanical loop the theory exhibits an elastic behavior and does not predict the 
tensile yield which is observed experimentally in Fig. 24. 

Figure 66 shows the thermomechanical loop predicted by Krieg., Swearengen and 
Rohde’s theory. The theory predicts a relatively small amount of stress relaxation 
at the maximum temperature hold period and the loop stabilizes in one cycle. In 
the tensile portion of the thermomechanical loop the theory exhibits tensile yield- 
ing with a slight stress dip at the top of the tensile portion of the loop. 

By changing the material constants by a small amount in the functional theory 
and in Krieg, Swearengen and Rohde’s theory, it is possible to produce pronounced 
stress dips at the top of the tensile portion of the thermomechanical loop similar 
to that predicted in Fig. 69. This predicted dip in the stress response at the max- 
imum tensile position in the thermomechanical hysteresis loop may be due to two 
factors. First, the strain rate at the maximum tensile stress is very small and the 
material may be relaxing. Experimentally the material exhibits a relatively rate- 
insensitive response at the temperature corresponding to the maximum tensile stress, 
so that experimental relaxation should not be, and is not, observed. However, the 
theoretical models cannot reproduce a rate-insensitive response, so that the 
theoretical models respond to a small strain rate by producing a stress relaxation. 

The amount and rate of relaxation depend on how far the stress at the maximum tensile 
position exceeds the equilibrium stress. Secondly, the theoretical models predict 
that during elastic deformation the equilibrium stress remains constant. During 
the initial compression phase of the cycle, the equilibrium stress is negative. At 
the maximum temperature of 954°C (lybO'^F) the hysteresis response starts to elastic- 
ally unload. As the stress increases elastically into the tensile region of the 
hysteresis loop, the equilibrium stress remains constant at its maximum negative 
value, since the equilibrium stress state variable changes only with inelastic deforma- 
tion. However, if the equilibrium stress increased due to changes in temperature 
as the stress increases elastically, the difference between the actual stress 
and the equilibrium stress 9 ,^^ would remain relatively small and the inelastic 
strain rate expression , which depends on the difference between and > 
would provide a gradual growth in the inelastic strain, similar to that observed 
in an isothermal hysteresis loop. But if ’’sticks*’ at its maximum negative value 
during an elastic excursion in which the stress is increasing, the inelastic 
strain will grow very rapidly. This rapid growth in inelastic strain will also 
produce a stress dip in the hysteresis loop. 



The equation governing the growth of the equilibrium stress in the 
functional theory may be written in the form 




On differentiation with respect to time, t, the following relation is obtained: 


iljj = (n,+n2)Cjj + 


dn, 

d® 




(6.5) 


If the material constant n£ (0(t)) appears inside the integral in the form U 2 ( 0 ( 0 ) 
the terms involving 0 in the preceding expression do not appear. However, with the 
©terms, the equilibrium stress can grow with temperature during an elastic excur- 
sion (in which c^^ = 0 and G = 0) according to the relationship 




dn, 

d@ 


® 


Cj| + 


I ^^2 

h2 d® 


® (fin 


0 ^ 
■fl jj — ri| Cjj ) 


( 6 . 6 ) 


In the FORTRAN subroutine HYPELA the functional theory has been implemented 
with the equilibrium stress having the form 

o "IgW) — G{£)1 <^ii 

n|(®($))“^dC +/ n2(@(C))e * dC . (6.7) 

so that the 0 terms will not be present in the corresponding differential equation. 
This form of the equilibrium stress equation prohibits the equilibrium stress from 
changing during elastic excursions and was used in the computation of the thermo- 
mechanical loop in Fig. 64. It is thought that the expression in (6.4), which has 
been written as part of the functional theory in Appendix 1, will give a better 
representation of the thermomechanical hysteresis response. Future investigations 
will implement Eq. (6.4). 

Stress dip anomalies due to relaxation at the maximum tensile position in the 
thermomechanical hysteresis loop can also be avoided if the instantaneous response 
of the material is assumed to be inelastic. At present the instantaneous response 
is elastic. If the instantaneous response is given by a^. = and the equilibrium 

response by , then by adjusting the growth laws for and so that 

at any time their magnitudes in the inelastic region are close to one another 
(5- • « . . ) , the response can be made rate- insensitive at lower temperatures and 

the^ relaxation in the tensile poriton of the thermomechanical hysteresis loop can be 
avoided . 
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Based on the favorable comparison of the experimental and theoretical hysteresis 
loops, together with good short time creep predictions, at temperatures ranging from 
800°F (427 °C) to 1800°F (982'*C), the functional theory was chosen to perform the 
cumbustor liner rig analysis. 


6.5 Combustor Liner Rig Analysis 

A three dimensional inelastic finite element analysis of the induction rig 
specimen described in Section 4.4 was undertaken in Ref. [8] under contract NAS3-21836 
with the NASA— Lewis Research Center . This analysis was conducted with constitutive 
models currently available in the MARC program. The analysis has been repeated in 
the present contract with the functional viscoplastic theory. 

A three dimensional finite element model of the cylindrical louvered induction 
rig specimen is shown in Fig. 67. The model contains an included angle of 0.577° 
which represents the angular distance between the radial symmetry planes of two 
adjacent cooling holes. This small model size was selected because of the circum- 
ferentially uniform temperatures in the louver lip region and to minimize computer 
run time. 

As shown, the model consists of thirty (30) finite elements. Twenty-one (21) 
of the elements are the twenty (20) node isoparametric element with reduced inte- 
gration using eight (8) Gauss integration points (MARC element #57). Since it is 
known that some reduced integration point elements are unreliable if distorted shapes 
are used, the full twenty seven (27) integration point element (MARC element #21) 
was used around the cooling holes. Selection of the reduced integration point ele- 
ments was based on initial estimates of the savings in computer run time over the 
full integration point element, and to minimize the required spatial temperature 
interpolation between the heat transfer and structural analysis meshes. 

The effect of the complete shell structure was simulated by the appropriate 
boundary conditions. Along the radial plane AC and BD in Fig. 67, only radial dis- 
placements were permitted. This was accomplished by use of the TRANSFORMATION option 
in MARC to transform the global degrees of freedom in these planes to the local 
coordinate systems. Along the planes AB and CD the effect of the fore and aft 
louvers was simulated by requiring that the radial displacements, u^, of nodal points 
on each plane were related by the equation 




lAflL 

^CD 


( 6 - 8 ) 


where ®-AB is the original radial coordinate on AB and is the original radial 
coordinate on CD. An additional condition equated the axial slopes along these 
planes. These conditions were prescribed by means of the user subroutine UFORMS. 
The fact that a series of louvers may be represented by this technique had been 
demonstrated in previous elastic shell of revolution analyses. 


54 



Numerical accuracy problems associated with the small included angle were 
investigated with a simple two element model having an included angle equal to 
that of the louver model • An initial analysis with this model produced stresses 
on the order of 69 MPa (10000 psi) when run in an isothermal condition in which 
the stresses should have been zero. In this model the radial displacements, due 
to the free thermal growth, were serveral orders of magnitude larger than the ^ 
circumferential displacements. Circumferential (hoop) stresses were produced as a 
result of numerical round-off during the single precision solution on the IBM 370/ 
3033 computer. Using a double precision version of MARC significantly reduced the 
round-off error to produce an acceptable isothermal solution, 

"The thermal increment file generated in Ref, [8] was used to drive the struc- 
tural analysis, A small mechanical pressure load was applied initially followed by 
sixty-one (61) thermal load increments per hysteresis cycle. Two hysteresis cycles 
were completed in the analysis. The transient temperature at the integration point 
closest to the end of the louver lip (the fatigue critical location) is shown in 
Fig, 68, Figure 69 depicts the hoop stress vs, hoop mechanical strain hysteresis 
loop at the corresponding integration point for two thermal loading cycles consist- 
ing of one hundred and twenty-one (121) MARC thermal load increments. Thirty (30) 
subincrements per MARC increment were used to ensure an accurate description of the 
viscoplastic behavior in the subroutine HYPELA. With the new integration algorithm 
the two thermal loading cycles could have been traversed with much less than one 
hundred and twenty-one (121) MARC thermal load increments; this number of increments 
was used, however, because of the availability of the thermal load increment file 
from the analysis conducted in Ref. [8]. 

The thermomechanical response shown in Fig, 69 exhibits an apparent ratchetting 
in the negative strain direction. This ratchetting is believed to be due to the 
elastic behavior observed in the low temperature tensile portion of the thermome- 
chanical hysteresis loop, A stress dip phenomenon is also observed at the maximum 
tensile stress obtained in the loop. Neither the excessive elastic behavior, nor 
the stress dip phenomenon, is observed in the experimental thermomechanical response 
shown in Fig. 24. Note, however, that no ratchetting is observed in the experimental 
thermomechanical response in Fig. 24 because the maximum mechanical strain amplitude 
is kept constant at -0.45 percent. 

The ratchetting of the thermomechanical hysteresis loop with its associated 
tensile elastic response and stress dip phenomenon, coupled with the generally 
poor predictions of the thermomechanical response produced by the viscoplastic 
constitutive theories, points to a need for further research on the thermomechanical 
behavior of Hastelloy-X and other high temperature structural materials. In general, 
the isothermal predictions of the theories with respect to strain rate, creep and 
relaxation behavior, is good. Further analysis, along the lines suggested in 
Section 6.4, is required to improve the thermomechanical predictions. 
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7.0 DISCUSSION OF RESULTS AND CONCLUSIONS 


The unified viscoplastic constitutive theories of Walker (Appendix 1) , 

Miller (Appendix 4), and Krieg, Swearengen & Rohde (Appendix 8), were installed 
for use in the MARC nonlinear finite element program by means of an initial stress 
method. An experimental program on Hastelloy-X material supplied the material 
constants for each of the theories at temperatures varying between 800°F (427 °C) 
and 1800°F (982®C) . General conclusions and observations concerning the work per- 
formed in this contract are: 

L. The. unified viscoplastic theories are represented by systems of '*stiff** 

differential equations. In order to integrate such theories over the finite 
element load increment, it is necessary to split the finite element increment 
into a number of subincrements and integrate the constitutive theories over 
the subincrements. By using a large number of sub increments the constitutive 
theories can be integrated accurately over arbitrarily large finite element 
load increments. The size of the finite element load increments is limited 
only by the fact that proportional loading (strain increments are proportional 
during the finite element load increment) is assumed to occur over the incre- 
ment. Consequently, the finite element load increments must be small enough 
to ensure that the strain history at any point in the structure can be repre- 
sented by increments of proportional loading. The strain history, at any point 
in the structure, which in general is nonproportional, is then represented by a 
series of proportional strain increments in which the constant of proportionality 
changes from increment to increment. 

The solution time required for a complete load increment in a finite element 
program is linearly related to the number of sub increments used in the inte- 
gration of the constitutive equations. A reduction in the number of subincre- 
ments can be effected if the viscoplastic constitutive theories are written in 
/an integral form and the theories are integrated over the subincrements by 
means of a suitable recursive algorithm. By means of the recursive algorithm 
developed in this contract the number of sub increments can be reduced to one. This 
corresponds to integrating the constitutive equations directly over the finite 
element load increment. Because the recursive algorithm is stable and accurate, 
even when the integration increment is large, the computer time required to 
integrate these "stiff” viscoplastic formulations can be substantially reduced. 

2. The viscoplastic constitutive theories examined give adequate accuracy in 

predicting the rate-dependent response of Hastelloy-X under uniaxial isothermal 
loading conditions. It appears that tabulating the material constants as a 
function of temperature provides improved accuracy over that obtained by assum- 
ing that strain rate effects and temperature effects are coupled by the Arrhenius 
type of exponential relation used in Miller's theory. Care must be taken, when 
tabulating the material constants as a function of temperature, to ensure that 
a reasonable asymptotic stress-strain response is obtained. Some of the material 
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constants, such as the inelastic strain rate exponent, n, control the magnitude 
of the computed stress at large strain values. The best way to obtain an inter- 
polated value of n is to interpolate the asymptotic stress-strain response and 
then to determine the value of n from the interpolated response. In tliis con- 
tract the reciprocal, n“^, is interpolated in the temperature tables, since this 
gives an adequate representation of the asymptotic stress-strain response a^ a 
given temperature, 

3, Relatively poor accuracy is obtained when the viscoplastic theories are used to 
predict the uniaxial thermomechanical behavior of Hastelloy-X, It appears that 
improved accuracy can be obtained by assuming that the state variables can change 

'' With temperature during "elastic” excursions and by assuming that the instantan- 
eous material response is inelastic, rather than elastic, 

4, The computer time required for the three dimensional structural analysis of the 
combustor liner induction rig test specimen was about the same as that required 
by the conventional yield surface plasticity theories available within the MARC 
program. This structural analysis was performed by integrating the functional 
theory over thirty (30) subincrements for every MARC load increment. Since the 
number of subincrements can be decreased by a factor of thirty and the number of 
MARC load increments can also be substantially reduced compared with the number 
of increments required to preserve accuracy and stability with the conventional 
yield surface plasticity theories, the new integration algorithm will cut the 
cost of nonlinear finite element analyses substantially over the cost required 
for conventional plasticity analyses. 
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8.0 RECOMMENDATIONS 


This contract has provided an enhanced capability for analyzing structural 
components which operate under cyclic elevated temperature conditions. Further 
work is necessary in order to: ' 

1. Improve the uniaxial thermomechanical response of the constitutive theories. 

Improvements can be effected by: (a) allowing the state variables to change 

with temperature during "elastic” excursions in a thermomechanical simulation; 
and (b) by allowing the instantaneous response, = “ij ’ inelastic by 

specifying suitable growth laws for the variable . 

2. Verify, and modify if necessary, the constitutive response under multiaxial 
loading conditions. 

3. Apply the constitutive theories to other materials, both isotropic and aniso- 
tropic. 
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Figure 1. Typical Combustor Construction Showing Predicted Temperature Distribution 
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Figure 3. Current Structural Technology is Limited by Inaccuracies in the 
Viscoplastic Constitutive Equation and Life Prediction Rules 


81 - 9 - 103-2 


64 

















• CREEP RATE DEPENDS 0N<7 — Q 

• a—Q, DECREASES DURING CREEP 


Figure 6. Creep and Relaxation Rates Depend on the “Overstress.” Different 
Behavior is Observed on the Loading and Unloading Branches of 
the Hysteresis Loop. 
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Figure 7. Comparison of Creep Response of Hasteiioy-X at 871 ®C (1600®F) on the Loading and 
Unioading Branches of a Steady State Hysteresis Loop Executed at a Constant Strain 
Rate of ±1.4 X 10‘^ sec'1 with a Strain Ampiitude of ±0.4% 
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Figure 8. Negative Stress Relaxation Response of Hastelloy-X at 871 ®C (1600®F) 
Initiated from a Steady State Hysteresis Loop Executed at a Constant 
Strain Rate of ±1.36 x 10’^ sec'^ with a Strain Amplitude of ±0.4% 
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Figure 9. Characteristic Tri-Llnear Stress-Strain Behavior Exhibited by those Theories 
(Miller; Hart; Krieg, Swearengen and Rohde) In which the Equilibrium 
Stress Grows Linearly with Inelastic Strain and Reaches a Saturation 
Limit in which Hardening and Static Thermal Recovery are Balanced 
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Figure 10. Strain Controlled (SC) Specimen 
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Figure 11. Experimental Steady State Hysteresis Loops for Hasteiloy-X Showing 
Inverse Strain Rate Sensitivity at 800 °F (427 **C) 
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Figure 15. Experimental Hysteresis Loops for Hasteiloy-X at 871 (1600°F) 
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Figure 16. Experimental Hysteresis Loops for Hastelloy-X at 
1800®F at ±0.6% Strain Amplitude 
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Figure 18. Theoretical and Experimental Creep Behavior for Hastelloy-X at 982X 
(1800°F) Starting from Various Points on the Tensile and Compressive 
Loading Branches of a Steady State Hysteresis Loop Executed at a 
Strain Rate of 3.66 x 10'^ sec‘^ at ±0.6% Strain Amplitude 
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Figure 19. Standard Tensile and Creep Sheet Specimen (Units in cm) 
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Figure 20. Theoretical and Experimental Relaxation Behavior for Hastelloy-X at 
982 (1800*’F) Starting at a Zero Strain Value on the Compressive 
Branch of a Steady State Hysteresis Loop Executed at a Strain Rate 
of 3.70 X 10'^ sec'l at ±0.6% Strain Amplitude 
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Figure 21. Theoretical and Experimental Relaxation Behavior for Hastelioy-X at 
982 X (1800**F) Starting at a Zero Strain Value on the Compressive 
Branch of a Steady State Hysteresis Loop Executed at a Strain Rate 
of 3.66 X 10'^ sec'* at ±0.6% Strain Ampiitude 
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Figure 22. Theoretical and Experimental Relaxation Behavior for Hasteiioy-X at 
982^C (1800°F) Starting at a Zero Strain Vaiue on the Compressive 
Branch of a Steady State Hysteresis Loop Executed at a Strain 
Rate of 3.87 x 10"^ sec'"* at ±0.6% Strain Ampiitude 
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Figure 23. Tubular Specimen for Thermomechanical Testing 
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Figure 25. Thermomechanical Temperature History 
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Figure 26. Thermomechanical Strain History 
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bustor Liner induction Heating Test Rig 
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Figure 29. Assembled Liner Test Specimen and Cover (a) and 5-Louver Test Section (b) 


Reproduced from 
best available copy. 


90 


81 - 9 - 103-15 



TEMPERATURE 


871 

( 1600 ) 

816 

( 1500 ) 

760 

( 1400 ) 

704 

( 1300 ) 

649 

( 1200 ) 

593 

( 1100 ) 

538 

( 1000 ) 

482 



CYCLI 


1 



Figure 31. Flow Chart of MARC Iteration Procedure 
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Figure 33. Theoretical Hysteresis Loops for Hastelloy-X Predicted with Linearly 
Interpolated Material Constants 





Figure 34. Axisymmetric Finite Element Model Used in the Thermomechanical History 

Simulation Models One Quarter of the Problem Using Roller Boundary Conditions 
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Figure 35. For Large and Small Work Hardening Rates the Stress-Strain Curves are 
Relatively “Square”. For Intermediate Work Hardening Rates the 
Stress-Strain Curves are Rounded 
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Figure 38. Theoretical Hysteresis Loops at 649 ®C (1200®F) 
Generated with Functional Theory 
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Figure 39. Theoretical Hysteresis Loops at 704 °C (1300°F) 
Generated with Functional Theory 
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Figure 41. Theoretical Hysteresis Loops at 816®C (1500“F) Generated 
with Functional Theory 
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Figure 42. Theoretical Hysteresis Loops at 871 °C (1600°F) Generated 
with Functional Theory 
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Figure 43. Theoretical Hysteresis Loops at 927 ®C (1700®F) Generated 
with Functional Theory 
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Figure 44. Theoretical Hysteresis Loops at 982 **C (ISOO'^F) Generated 
with Functional Theory 
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Figure 45. Theoretical Hysteresis Loops at 649 ®C (1200®F) Generated 
with Krieg, Swearengen & Rohde’s Theory 
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Figure 46. Theoretical Hysteresis Loops at 704 °C (1300°F) Generated 
with Krieg, Swearengen & Rohde’s Theory 
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Figure 47. Theoretical Hysteresis Loops at 760 ®C (1400®F) Generated 
with Krieg, Swearengen & Rohde’s Theory 
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Figure 48. Theoretical Hysteresis Loops at 816°C (ISOO'^F) Generated 
with Krieg, Swearengen & Rohde’s Theory 
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Figure 50. Theoretical Hysteresis Loops at 927 °C (1700“F) Generated 
with Krieg, Swearengen & Rohde’s Theory 
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Figure 51. Theoretical Hysteresis Loops at 982 (1800°F) Generated 
with Krieg, Swearengen & Rohde’s Theory 


81 - 9 - 103-27 


112 



STRESS (ksi) 1 ksi = 6.9 MPa 



Figure 52. Theoretical Hysteresis Loops at 649 °C (1200‘’F) Generated 
with Miller’s Theory 
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Figure 54. Theoretical Hysteresis Loops at 760 ®C (1400“F) Generated 
with Miller’s Theory 



115 


81 - 2 - 112-4 







STRESS (ksi) 1 ksi = 6.9 MPa 



STRAIN (7o) 


Figure 57. Theoretical Hysteresis Loops at 927 ®C (1700®F) 
Generated with Miller’s Theory 
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Figure 59. Comparison of Experimental and Functional Theory Predictions for 
Hastelloy-X at 871 °C (1600'»F) 
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Figure 62. Hysteresis Loop Predicted by Functionai Theory at 982 (1800°F) for Strain Rate of 
1.25 X 10'^ sec*'*. Three MARC Increments were used to Load from Zero Strain to 
0.6% Strain with Thirty Subincrements per MARC increment 
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Figure 64. Thermomechanical Loop Predicted by Functional Theory 
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TEMPERATURE (°F) 


1800 



Figure 68. Transient Temperature Response at Integration 
Point Closest to the Louver Lip 
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Figure 70. Elastic-Perfectly Plastic Hysteresis Loop Approximation 
for Hastelloy X at 982^*0 (1800“F) 
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Figure 71. Temperature Rise in Hastelloy-X Specimen at 982*’C 
(1800*’F) Due to Plastic Deformation 
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Appendix 1. Walker’s Theory (Integral Form) 




dei 




'I _ ^ 2 <^_kk _ ^ j ^ _ ( 1 ) 


njj{t) = nij(t) + n,(®(t))Cij{t) + n2(©(t))^'e“i®‘^’"®‘^’( d^, 


( 2 ) 


K(t) = K|(@(t))-K2{®(t))e""7<®(t))R(t)^ 


(3) 


LO 

U> 


= J ^SjjX(®(t))-J^ +2^(®(t)) ^ Sjj a(®(C)) ^3\(®(C))+ 2;x(®(C))j d£, 


c • • (t) = 

^ \^IJ 
0 




0 Q 

::Xl(®(t))+3n(®(t)) T-r , 




(4) 


(5) 


/•♦ 3/x(®(C)) /aR \i-i/n(@(f)) 


G(t) 


=_^'{(n3(0 


m ( 0 (g)) -I 

(®(C)) + n^(®(C))e“"^^®'^”^‘^’ j-||- + n6(®(0) (O ^ 


' }dC, 


( 6 ) 


(7) 


R(t) = 


t 5 Cjj dc ij 



'oV ^ dc ac 


dC 


( 8 ) 


Material constants: X, /x, II ,n,m ,n,, 02,03,04, ng.ng.ny.K, ,K2 depeod on temperature. 



Appendix 2. Walker’s Theory (Differential Form) 



( 1 ) 


hji = (n| + n2)C::+C 


(3n 


'J d@ 


® -(iijj -iijj-n 


i^ij) 


j-IHl 

nz a® 


® 


( 2 ) 


K = K 


K o 6 


— n-^ R 


(3) 


Cij = (SjjXe,^,^ + 2^ejj-ai. - Sjj(3X + 2^) a® )/2^ , 


(4) 


6 = 


(03 + 04 e 



XI-. D.-- 
'J 'J 



(5) 


R = 



CijCjj , 


( 6 ) 


I » 

■ “^ij 3 °ij 

0 

Material constants: X,/i ,n,m ,n, , 02 , , 04 , Og, Og , , K , , K2 depend on temperature. 



Appendix 3. Chaboche’s Theory 



ai, = cf(R)[aCi,-Rflii-b (-f + 

af(R) 

if GijOii)+giR). 

Cij = ^SjjX€kk + 2;xeij-cr,.-Sij(3X + 2^)a®)/2/i , 
f(R) = i + (l-i)e 


g(R) =h(i-e , 



^■kk' 


•^ij) 

i)(iSjpi2j|) 


J_ + _L 
a d® f 


Material constants; X,/i,K,k,n,m,a,b,c,i,h,;9,/ depend on temperature. 



Appendix 4. Miller’s Theory (Differential Form) 



k=HgR ^Cg + ^-f ilijiijj' - -^K^j-HgCgB^'IsinhlAgK^)}" , 


Cjp (Sjj \k + 2^€jj -6-jj -Sjj(3X+2^)a@)/2/z 

H-* 

00 

ON 

^ij ~ ^ij^’kk ’ 

o'- exp ■! !■ for T > 0.6Tfin , 

e'-exp for T< 0.61^ 


(1) 

( 2 ) 

(3) 

(4) 

(5) 

( 6 ) 

(7) 

(8) 


Material constants: n ,H, , H2 , A| , A2 , C2-iO*> k are independent of temperature. 

Material constants: X, ,Kq ^O' depend on temperature, where is the initial value of K, and T is the temperature in °K. 



Appendix 5. Miller’s Theory (Integral Form) 


o 








(1) 


LO 






-)G(t)-G(C)( acij 


ac 


dC. 


K(t) = Ko + H2 ^C2+y|-nj.(C)nij(C)‘ -A2K3(C)/A,je , 


Cij(t) = 


0(t) 


/ a€^. a€ij da-\ 

f' , . -,l I aR 

Jq k(C) \Be'(©(C)) ac 


a^ij ao-ji / \ a® \ 

_j — SijQ(®(C))(3x(®(e))+2Ai(@(C))l-^l dc , 


" ( 3 aR 

ac dc, 


t {sinh(A,y-|-nij(emij(o )}" 

G(t)=H, / B0'(®(C)) — di, 

° 

/.t , {sinh(A2K^(C))}" 

j(t)=H2C2 ^ Ba'(@(C)) ^ ■ dC , 


( 2 ) 


(3) 


(4) 


(5) 


( 6 ) 


(7) 



Appendix 6. Lee and Zaverl’s Theory (Isotropic Form) 



K = qY, 

CjjMSijX€|^l^ + 2^€ij-ajj-Sij(3X+2^)a@)/2^ ’ 


Material constants: X ,/x ,Kp,YQ,n ,q ,C, g, h.a depend on temperature, where Kq and Yq are the initial values of K and Y. 
Material constants: B ,D,k are temperature independent. 



Appendix 7. Bodner, Partom and Stouffer’s Theory (Isotropic Form) 





ij(C) 

Cij(0! 




r* Af 

Z = m(z,-ZQ)e dC . 

Material constants: X,/x ,Dq .z,q,m ,n ,Z|,Zq depend on temperature 

The quantity Zy is not a tensor, but depends on the sign of the ij^*^ component of the inelastic strain rate. 



Appendix 8. Krieg, Swearengen and Rohde’s Theory (Differential Form) 


^ii 


\/|('l'Si) ^ij)(^^i) ^ij)* 









j 

pq^pq 


(• 


2 



K^A^R-AglK-Ko^ , 


R=./tCijCjj 



'^kk ■ 


Material constants; X,;j.,n ,A,,A2,A3 ,A4,A5,K(^ depend on temperature, where is the initial value of K. 



Appendix 9. Krieg, Swearengen and Rohde’s Theory (Integral Form) 








^kk 2 


K(t)=K(,(®(t)l*^A3(e(0)e 


dC, 


de 


Cii(t)=j ( SiiX(®{C))-^ +2/x(@(e)) 


deij < 5 o-j 






ae dc 


-8ija(®(0) (3X(®(e)) +2^(®(C)))-f|-j dC, 


Q(t)= f 


' 3/x(®(C)) / aR \'-'/n(®(0) 

'o K(0 lac) 


dC, 


/? '( A3(@(C))-|-nii(C)n::{C) \ . 

J(t)=j^ A4(®(0)(K(C)-Ko(®{f)))"‘®‘^”"' dC , 


( 1 ) 


(2) 


(3) 


(4) 


(5) 


( 6 ) 


(7) 


R(t)- /' /2 acjj acij~ * 

•'ov 3 ae ae 


( 8 ) 
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Appendix 10. Cernocky and Krempl’s Theory (Integral Form) 




( 1 ) 


'1 


( cosh b(@(t)) +a(© (t) ) ^ M ..(t)'i': :(t) 

V ^ 


( 2 ) 


xPi.(t)-- 


SjjX(®(t))€|^k(t) + 2^(®(t))€jj(t)-Sjja(®(t))|3X(©(t))+2^(®(t))| 


/E(®(t)), 


(3) 



j(t)=K(®(t)) exp 


S(® (t))exp 


C (®(t) )y(cTij(t)-ilij(t))(cr.j(t)-llij(t))' 


fj. {'5\+ 2 fj.) 

Material “constants”: X,/i.,Epa-,b,K,S)^ depend on temperature, E = Young s modulus = \ + ij_ 

Material “constants”: E^,a,b,K,S are updated under cyclic loading conditions. 

In this particular formulation it is assumed that Poisson’s ratio is constant. This assumption can easily be relaxed. 


(4) 


(5) 
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Appendix 11. Hart’s Theory 




Cij =(Sjj + 2^€jj-CTj. - S|j(3X+2/A)a©)/2^ , 

Material constants: X,;x, at K,n,f ,o-Q,nn , /3,C ,/,S depend on temperature where a* is the initial value of cr* 
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Appendix 12. Valanis’ Theory (Isothermal Form) 


ci.(t): 


Sii 


|x+ G [z(t)-z(C)] 


I <36 






Sii 


kk 


dC 


( 1 ) 


r T -a.[z(t)-z(C)l -02fz(t)-z(^)| 

G[z(t)-z(0] = G,e ^Gge J 


( 2 ) 


z(t) = -^i/^(|-/3R(t)), 


(3) 


it ' (it) • 


(4) 


ae 4 <3cij acij 


ac V 3 ac ac 


a function of which varies [25] according to the type of deformation (creep, relaxation), 

V / ac 


(5) 


( 6 ) 


aCii 


de 


= S„X^+2,.-^-k-^ ,(0<1«I), 


aeii a 


o-ij 


(7) 


-^ir ac ■ ac “ ac 

Jm.) 

Material constants; X,yu., G, ,62 ,Q|, “2 >/3 >*'• Further constants are needed to define T\^a^/ 



APPENDIX 13 


A question naturally arises as to how much heat is generated in uniaxial ^ 
Hastelloy-X specimens when they are subjected to cyclic, triangular s^train his- 
tories, at strain amplitudes of about ^ 1 percent. The specimens are cycled at 
a given temperature which is monitored by thermocouples attached to the grips at 
the ends of the specimen. If heat is generated by the dissipation of inelastic 
work during cyclic loading of the specimen, the temperature of the specimen will 
increase above that of the oven temperature until a steady state is reached. 

Under steady state conditions, the rate at which heat is produced by inelastic 
working of the cylindrical specimen is equal to the rate at which this heat is 
conducted through the grips at the ends of the specimen and through radiation and 
conduction at the cylindrical surface of the specimen. The following calculations 
show that at 982°C (1800°F) , for a strain rate magnitude of 3.9 x 10 ^ sec , the 
steady state temperature in the central part of a Hastelloy-X cylindrical button 
head bar specimen exceeds the oven temperature by about 0.8°C (1.4®F) due to in- 
elastic working. At room temperature, the temperature rise in the specimen does^ 
not exceed 22°C (40°F) at a strain rate of 3.9 x 10 ^ sec > whilst at 3.9 x 10 
sec ^ it does not exceed 2°C (4°F). Hence, the effects of inelastic dissipation 
of work into heat can be neglected in the data reduction of the Hastelloy-X tests. 


The coupled heat conduction equation for a cylindrical rod under uniaxial 
loading may be written in the form 


ps 


)f2 '\dr^ 


I a® 

' dr az 


@ \ 

-2- j+D , 


( 1 ) 


where p is the density, s is the specific heat, 0 is the temperature in the speci- 
men, and D is the rate of dissipation of work into heat due to inelastic working. 

In formulating Eq. (1), small heat source terms, such as that arising from the ex- 
pansion and compression of the material, have been neglected in comparison with 
the inelastic dissipation. 


Under cyclic loading, about 90 percent of the work in typical metals is dissi- 
pated as heat whilst the remainder is used to produce microstructural changes 
(evidenced by working hardening) . Since Hastelloy-X does not work harden at 
elevated temperature, one may assume that all of the work is dissipated as heat and 
take D in the form D = a c where c is the inelastic strain rate. The average dis^ 
sipation rate over a hysteresis cycle is then obtained in the form 



where T is the time for one cycle. 


( 2 ) 
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The value of D may be found by noting that the stress-strain hysteresis loop 
for Hastelloy-X at 982°C (1800®F) may be represented fairly well by an elastic- 
perfectly plastic loop. Hence, 

T 

I / _ orea of loop in fiq. 70 

T — — 


or 


D = 






where a is the maximum stress in the hysteresis loop, 
max 

At large strain values (^1%), the total strain is approximately equal to the 
plastic strain, e « c, so that the magnitude of the total strain rate is given by 
|e| « 4A/T. The average dissipation rate then becomes 




(4) 


Since a and |e|_are spatially constant throughout the rod, the average dissipation 
rate, given by D in Eq." (4), represents a constant source of heat. 

The steady state equations governing heat conduction in the rod are: 


a^@ 


d^@ 


. D 




k-|^ + h(0-©o)--O at na, 


@:®o at z =0 and z -a , 


(5) 

( 6 ) 
(7) 


where is the temperature of the end grips and the oven (as measured by thermo- 

couples on the grips) . 

Choosing 


®=®o + v(r,z)+<Jb(z), 


( 8 ) 
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Eqs. (5)-(7) become 





f + '^'t+0=0. 

(9) ^ 




dr^ r dr dz^ 





k +hv= h0atr = a, 

(10) 




\/--4> at 2=0 and z-f . 

(11) 

Equation 

(9) 

can 

be made homogeneous by choosing 







Choosing 

= 

0 at 

z = 0 and z = i then gives 





<^ = ^Q2(i-z), 

(12) 

and Eqs. 

(9)- 

-(11) 

become : 





<3^ 

dr^ ' 

(13) 




k~^~ + hv= ^hQz(i— z) at r-O) 

(14) 




v=0 at z =0 and z-i . 

(15) 


The solution to Eq. (13) which is finite at r *= 0 and satisfies the boundary condi- 
tion in Eq, (15) is 


? n Hz / n nr 

v = X On sin i To\ i 


where I^(x) is the Bessel function of order zero with imaginary argument given by 


Iq(x) = Jq( j x) ~ Y, 
n=0 



(16) 
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Satisfaction of boundary condition (Eq. (14)) requires that 



I,(x) = 


dlo(x) ^ 
d X n=0 



n+i 


n !(n+i) ! 


(17) 


The Fourier expansion of z(Z-z) on the interval (0,^) is given by 


00 

zU-z)- Z Bo 

n=i 


nnz 
sin g 


where 



z, 


or 


00 

z(i-z)=^ 

n=i 




{(-l)"‘^'+l} sin 


nnz 


(18) 


Comparing coefficients in Eqs. (17) and (18) gives 


2hQi^ 


1)0 + 1 

- 


n3 j 

knO T j 

f -^ol 

1"5“1 

l+hlol 

(^) 


(19) 


For n = 2,4,6,... the coefficients a^ in Eq. (19) are zero. Equations (8), (12), 
(16), and (19) then give the temperature in the rod in the form 




4hQf^ ^ 


. (2n+onz _ r(2n+i)nf 
sm — j IqL j 


nfo(2n + l)^ | M 2 n + i)n |~ iinyriQ ^ 


(2n-H)no 


■]} 


( 20 ) 


The maximum excess temperature occurs in the middle of the rod where r — 0 and 
z = £/2, and is given by 
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®-©o = A® - 


8 


lhQii_y I (-1)^ 

n3 (2n+l)^ |~ k(2n+i)n |~ (2n+nna' , • (2n+i)no j| 


From the relation Q = D/k = a lei, this becomes 

max 


_ 4ho-maxl gU^ 


8k 


kn 


kfl T 

~r^\ 





( 21 ) 


The thermal conductivity k for Hastelloy-X at 982°C (1800®F) is 


k = 195 


Btu. in 
ft?hr.°F 


195 X 778 X 12 
144 X 3600 


. 2 
in. lb 

2 

in. sec.°F 


3.51 


lb 

sec.^F 


The value of h is given by 


1 -lor Btu 135 X 778 X 12 in. lb ^ , lb 
h « 135 — = — = 2.43 

hr.ft.°F 144 X 3600 in.sec.°F in.sec.®F 

/ 

~3 ““1 

At a strain rate of \k\ = 3.866 x 10 sec (10 cycles per minute at - 0.6 
percent strain amplitude), <^^ax ” 19212 psi. The length of the specimen is £=1 in. 
and its radius is a = 0.15 in. 

Substitution of these values into Eq . (21) gives the excess temperature at the 
center of the rod as 


149 



From this result it can be seen that only the first term in the infinite series 
contributes significantly to the result. Figure 71 shows the computed steady state 
temperature excess, A0 , in various parts of the specimen obtained by keeping only 
the first term in the series expansion in Eq. (20). The radial distribution of ^ 
temperature in the bar is almost uniform and the axial distribution is the difference 
between a parabolic and sine function. 

At room temperature = 100,000 psi at 10 cycles per minute and k = 60 Btu 

in./ft^. hr. °F. The constant h is probably small, ~1 Btu/f t^ .hr . °F. Neglecting, 
the* series terms which account for the radial outflow of heat, one obtains A0 ~ 

45°F at 10 cycles/minute and A0 « 4.5°F at 1 cycle/minute at room temperature. 
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APPENDIX 14 


RECURSIVE INTEGRATION METHOD FOR FUNCTIONAL THEORY 


The uniaxial form of the function theory in Appendix 1 may be written, for 
isothermal deformations (0=0), as 

Jl(t) = S + n|C(t)+n2 

(5) 

( 6 ) 

(7) 



The decreasing exponential function in Eq. (1) implies that most of the con- 
tribution to the integral comes from the region about the upper limit, t, where 
Q(0 w Q(t). For general integrals of the form 




(3m(0 


dC, 


( 8 ) 
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it may be shown that the asymptotic expansion of the integral for large values of 
t is given by the relation 




[mQ + m) 


it 

(9) 


where dots denote differentiation with respect to time. 

The first term of Eq. (9) may be applied to the integral in Eq. (1) to give an 
asymptotic expression for the stress well into the plastic region at large values 
of strain and time in the form 


<T -Q. 


Ee-fl 

+' : 

Q 


( 10 ) 


where 


Q = 



dc 

E 

d€ 1 dcr 

\at j 

at 

■ K 

at E at 


l-i/n 


( 11 ) 


Since e >> o/E and Ee >> Cl in this large strain region, the asymptotic expression 
for the stress takes the form 


o- = n + Ke'^'’. 


( 12 ) 


This relation can also be derived by assuming that e >> a/E and Eg >> Q in 
the uniaxial differential equation appropriate to monotonic loading, viz., 

C=€-^=(^)" (13) 

In a finite element solution the strain increment Ae is supplied to the sub- 
routine at the beginning of each load increment, and it is assumed that the strain 
rate Ae/At during the load increment is constant. It is, therefore, constant during 
each sub increment in a multi-step technique. 
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The integrals in Eqs. (1) and (2) contain the expressions (E 3e/95 ~ and 

(3c/3C = 3e/35 - 3a/E3£;) . In an asymptotic expansion of these integrals (valid for 
large values of t) the strain rate 3e/35 is the dominant term in the two preceding 
expressions. Since the strain rate term is constant during the finite element load 
increment, it is therefore permissible, in an asymptotic expansion of the integrals, j 
to assume that the expressions (E 3e/3C ~ 3f2/30 and (3c/35) are constant during 
the integration. 

t 

Consider the value of the stress at time assuming that the values ' 

of the stress and equilibrium stress are known at time t^^. Then, from Eq. (1) 








= ( cTn - 6 


-{on + |-Qn[ 






(14) 


For small time increments the exponential term in the integral is approximately 
equal to one and the integrand (E 9e/35 “ 3^/35) tnay be assumed constant in the 
integration over the small time interval. It is also permissible to assume that 
(E 3e/3g - dQ/dg) is constant in the asymptotic integration appropriate to large 
time increments. Hence, for both large and small time increments, the integral 
may be approximated in the form: 



where the subscripts n and n4-l refer to the values of the variables at the times 
tn and 

The preceding integral may be written 





doit), 


(16) 


153 



and s inc e 


dQ(C) _E 3c 

dC ’ K 3| 


i“l/n 


may be assumed to be constant in the integral for both large and small time incre- 
ments, the integral on the right-hand side of Eq. (9) takes the form: 



Equation (15) may thus be integrated in the approximate form: 




€p)- 



(17) 


For large time increments the exponential term vanishes and the integral takes 
the asymptotic form (E Ae - A^2)/AQ = (Ee -^^)/0, in agreement with the result in 
Eq. (12). For small time increments the integral assumes the approximate value 
(E Ae ^ ^9.) . 

Substitution of Eq. (17) into Eq. (14) gives the recurrence relation 


°‘n+i-'^^n + l+(cTn“iin^e ^ + [tfen+i-^n )-(i^n+i “iin )] 


l_g“K+i Onf' 
Qn+l “Q n 


. (18) 


Similarly, the recurrence relation for the equilibrium stress integral in Eq. (13) 
may be written in the form: 


0 0 
•^n+i ' h| Cn + 1 + (fin~^~h|Cn) 6 


“|Gn+i“Gn| 


-e 


“|Gn+l“®n[ 


i-ngCCn+i Cn)l 


(19) 
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In the recurrence relations in Eqs. (18) and (19), the values of 
Gn+i at the end of the increment may be determined by the approximate relations 

- A - e: / . . E 

Qn+i Qn'^Q' K yAt j K 

and 

Gn+i ~ Gn " f *3 + 6 n + ij AR + ng|iip|^i 


m- 


At, 


( 20 ) 


Ac 


At 


i-i/n 


At 


(19) 


where AR and Ac are defined by the relations 


Cn+i~On = Ac=(€n+|-^n)~ 


(g-n+i-o-n) 

E 


( 21 ) 


and 


Rn+|-Rn = ARHAcI • 


( 22 ) 


Evidently, the values of and depend on | Ac | and are therefore 

implicitly defined in Eqs. (19) and (20), since | Ac ] depends on ®n+l 

depends on An iterative technique is therefore required to resolve these 

equations. This iterative technique is described in Appendix 15. 
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APPENDIX 15 


RECURSIVE INTEGRATION METHOD FOR FUNCTIONAL THEORY 
USED IN SUBROUTINE HYPELA 


The integration algorithm developed in Appendix 14 is now used to integrate 
the three dimensional equations of Appendix 1. The solution to the implicit re- 
cursive relations is obtained by means of a Newton-Raphson technique. 

In the following algorithm the symbol in front of a variable refers to 

the difference in the variable between the beginning and end of the increment. The 
time at the beginning of the increment is denoted by t and at the end of the incre- 
ment by t + At, 

Since the stress tensor for non-polar media is symmetric, the constitutive re- 
lations are written in vector form, so that = ^11 > ^2 ~ ^22 > ^3 ” ^33 > ^4 ” ^12’ 
*^5 ” ^23’ ^6 ~ ^13 similar definitions applying to and • 

In the following algorithm the symbol e^(t) refers to the mechanical strain, 
so that 


A€j(t) = A€|(t)-SiaA@,8i = 


1 if i < ^ 
0 if i >3 


where, Aej^(t) is the total strain increment. It is assumed that Qj^(t), e^(t) , c^(t), 
r^^(t), R(t) , Q(t) and G(t) are known at the beginning of the increment. Initially, 
at time t = 0, a»(t) = (t) = c • (t) = 0 and- R(t) = Q(t) = G(t) = 0, whilst ^2.(t) = 

To integrate the constitutive equations from time t to time t -h At, the follow- 
ing operations are required. 

Step #1. Compute the values of the material constants at the temperature 
0 (t + 1/2 At) corresponding to the middle of the sub increment . 

Step #2. Assume an initial guess for the inelastic strain increment Ac^. 

In a subincrement method this guess is taken equal to the devia- 
tor ic strain increment in the first subincrement of the current 
load step. In subsequent subincrements the guess is taken from 
the preceding sub increment of the current load step. 
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step #3. Compute the inelastic strain at the end of the subincrement from 
the relation 

Cj(t+At)=C|(t) + AC| 

Step #4. Set At) = as a first guess for the equilibrium 

stress at the end of the sub increment . 

Step #5. Compute the initial values of the equilibrium stress at the end of 
the subincrement from the relations 

R|(t + At)= |c^(t + At)+c|(t+At) + C 6 (t+At)}/Rn , 

fl, 2 (t+At) =-R + 2H{c|(t+At) + C 4 (t + At)+Cg(t+ At)}/iin , 

S3(t+At) = -ft + 2S|c|(t+At)+c|(t + At) + c|(t + At)}/Sii- , 

^ 4 ( 1 + At) = 2 n {c|(t+At)c 4 (t+At) + C 2 (t+A ;-)C 4 (t+At) + C 5 (t+At)C 6 (t+At)}/% , 

^g(t+At) = 2 a{c 4 (t+At) C 6 (t+At)+C 2 (t+At) C 5 (t+At) + C 3 (t+At) C 5 (t+At)|/^n , 
n 6 (t+At) = 28 . {c|(t+At)C 6 (t+At)+C 4 (t+At) C 5 (t+At)+C 3 (t+At)C 6 (t+At)}/S;^, 
where 

Sjj = -|-{cht+At)+c|(t+At)+c|(t + At) + 2 [cht+At) + c|(t+At)^-c|{t+At)]}• 

Step # 6 . Compute the cumulative inelastic strain increment from the relation 

AR = 


If /3i(ACi 


where 


•I if i < 3 


■2 if i >3 
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step #7. Let subscript m refer to the number of iterations. 

Set m = 0. 

$ 

step #8. Set m = m + 1. This is the start of the iteration loop. 

Step #9. Set = AR. Compute the cumulative inelastic strain at the 

end of the subincrement from the relation 

R(t+At) = R(t) + AR 

Step #10. Evaluate the drag stress at the end of the subincrement from the 
relation 


K(t+At) =K|-K2 e 


-nyR (t + At) 


Step #11. Evaluate the mth guess for AG from the relation 


. r -Dc R(t+- At)l 2 2 , . 1 

A CpD = 1^03 +ri4 e J AR+Og At /Sj j (t + At) 

^ i = i 


m-l 


Step #12. Compute the inelastic strain at the end of the sub increment from 
the rel'ation 


Cj(t+At)=Cj{t) + ACj 


Step #13. Compute the equilibrium stress at the end of the subincrement from 
the recursive relation 


fli(t+At , AGm)=Sj(t+At)+n|Cj(t+At) + (Aj(t)-Si(t)-n|Cj(t)) e +ACj(^^^-^ 

\ AGm 

Note that this is an implicit relation for the determination of ^i^(t + At) since 
AGjj^ depends on Sl^(t + At) if the recovery term involving n^ is present (cf. Step 
#11) . If n^ = 0 and there is no recovery, this relation for J^£(t + At) is explicit 
and no iteration is required for its determination. 

Step #15. Compute equilibrium stress increment from the relation 



Ailj = Aj(t+At, AGm)-^ii(t) • 
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If iig = 0, skip the Newton-Raphson iteration for + At) and 

go to step #20. 

Step #16. Compute the function f(AG^) where 

f (AGm) = AGm - [n3+n4 AR-ng At [f ^ /3- (t+At , AG^)]” 

L -J ^\z{ 

Note that if AG^, the guess for AG computed in Step #11, is exact, then f (AG^) = 0. 

In general f(AG^) ^ 0 and a Newton-Raphson iteration is required to improve the value 
of AG^, (Note that the exact AG^ satisfies both Step #13 and the relation f(AG^) = 
0.) 


Step #17. Compute the derivative 




[f /3j^j(t + At,AG 


m 


(3H j(t +At, AGp^) 
’ a(AG„) 


where 


— AOpn —AOr 


axij(t+At, AG^) 0 -AGm i Aa 

^ ^ =-(ni(t)-n..(t)-n,c,(t))e ^^'^+n 2 Ac./ 


d(AG 


m' 


i/ j \ I / Ub j \ 


J 


'J 


(AGm) 


by differentiation of the relation in Step #13. 


Step #18. Refine AG^ with a Newton-Raphson iteration 


f(AGm) 


AGm+rAGm 


3(AGppj) 

Step #19. Compute refined value for the equilibrium stress at the end of the 
subincrement from the relation 


ii.(t+At,AGrn+!) =a.(t+At)+n,Ci(t+At)+(fij(t)-iij(t)“n| C;(t))e ^ l±ie 


:Ag 


m-hl 


j\l/ I l| j \ 


ag 


m+l 
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Compute equilibrium stress increment from the relation 


f 


Ailj=Uj(t + At, AGm + i )-iij(t) . 

Step #20. Compute the mth guess for AQ from the relation 


3^ 

AQm= K(t + At) 



Note that AQ^ depends on aj_(t + At) through its dependence on AR (cf. Steps //’s 
22 and 24)- 

Step #21. Compute the stress at the end of the .sub increment from the recursive 
relation 


cjj ( t+ At , AQp^ ) - 2 


i(t + At) + Qj|^ o-j(t) + ^X+-|-^Z Afj'j 


•^o-j fij(t)- a j cTj(t)^e +{2 


-Aq 


-jcTjWje +(2 mA€;-=, I 


m 


where Qj =: 


■\ ifi < 3 


'0 if i>3 . 


This is an implicit relation for Oj^(t + At, ^Q^n) > since AQ^ depends on a^(t + At), 
which must be resolved by Newton-Raphson iteration. 

Step #22. Compute the inelastic strain increment from the relation 


AC}(AQrr^) =|ajX X A6j + 2/x A€| -cTj(t + At, AQrn)+o-j(t)|/2/x . 


Step #23. Compute inelastic strain at the end of the sub increment from the 
relation 


C|(t+At)=C j(t)+ACj(AQm), 
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step #24. Compute the cumulative inelastic strain increment 


AR(AQm) 



/3j ACj(AQm) 


Step #24. Compute the function F(AQ^) where 


F(AQm)= AQpf^ 


3^ 

K(t + At) 


AR(AQm) 

At 


i-l/n 

At . 


Note that was evaluated in Step #20 and is used to evaluate AR (AQjj^) in Step 

#24. If the exact value of AQ^^^ is known, then F(AQj^) = 0. In general, F(AQjjj) 7 ^ 0 
and AQ^ must be refined by a Newton-Raphson iteration. (Note that the exact AQ^ 
satisfies both Step #21 and the relation “ 0.) 

Step #26. Compute the derivative 


aF(AQm) 

-1- 1 

\ |\[ar(AQJ 

-i/n 

a(AR(AQm)) 

a(AQm) 

■' K(t+At) \ 

^ n ;[ At J 


d(AQpr|) 


where 


a(AR(AQr^)) 

^(AOm^ 


d 

<5( AQ|x|) 



The derivative 8 (Ac ^^(AQ^) ) /9 (AQ^) is obtained from differentiation of the 
relation for ACj_(AQ^) in Step #22 in the form 

a(ACi(AQrri)) I aa|(t+At,AQ^) 


a(AQp^) 2^ a(AQ|Tn) 
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Step #27. Refine AQ^ with a Newton-Raphson iteration 


^^m+i ’ ^Qm 


F(AQm) 

dF(AQm) 

a(AQm) 


Step #28. Compute refined value of stress at the end of the sub increment from 
the relation 


(T./t + At,AQm+|) =X^i(t+At)+Qi|i^a.(t)+|x + -|-^^^ A6jj 


+ l(Jj (t) ^ ilj (t)"“Q ; 




+(2;.A6i-ai47-I Aep^Ailij 


i-e 


-Ao 


m + i 


AQ 


m- 


Step #29. Compute refined value of inelastic strain increment from the relation 


ACj(AQrp+|)- '^QjX X ACj +2^ A^j — cr. (t+At , AQpp+i ) + o’ /2^ 


Step #30. Compute refined value of cumulative inelastic strain increment 


AR Z 3 [AC|(AQm+|)] 


Step #31. If increment is elastic after two Newton-Raphson iterations, exit 
from the iteration loop: 


If AR< 


I 

iOO 



m=2 , go to Step #35 . 


Step #32. If m = maximum number of iterations allowed, go to Step #35. 
Step #33. If m < maximum number of iterations allowed, go to Step #8. 
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step #34. If 


> 


— , go to step #8. 


AR - AR 
m 

AR 
m 

Step #35. The iterative loop is now complete. The state variables at the end 
of the sub increment are now given by the values of the quantities 
0(t + At), R(t + At), c^(t + At), Q^(t + At), e^^Ct + At), a^(t + At) 
at the updated time increment t -t- At. 

Step #36. Go to Step #1 for next time increment. 
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GO TO 903 
816 CONTINUE 

DO 821 J=l ,6 
6 ( J ) = SUMIN ( J ) 



Appendix 17. Fortran Listing of Functional Theory in Subroutine HYPELA (Continued) 




621 CONTINUE 

DO 819 J-1 ,6 
DO 819 Kzl ,6 
D (J,K)=0. 

819 CONTINUE 

DO 820 J:: 1 , 3 
DO 820 K-1 ,3 
A LPHA::G . 

IF( J.EO .K ) ALPHA-1 . 

D (J,K)-C5+ALPHA»C3 
32Q CONTINUE 
D (4,4 )-C4 
D (5,5):::C4 
D (6,6)=C4 
903 CONTINUE 

OMPUTL MARC STRESS INCREMENTS 
DO .822 J-1 ,NGENS 
SUM-Q. 

DO 823 Krl,NGENS 
SUMzSUM+D( J,K )^^DE (K ) 

823 CONTINUE 


DS( J)=SUM + G( J) 

822 CONTINUE , 

STATE VARIABLE INCREMENTS INTO STORAGE ARRAYS FOR NEXT MARC INCREMENT 
DTEMP ( 3 )rRB-TEMP ( 3 ) 

DTEMP( 16) ro. 

T£MP(16)-AK1-AK2>:‘EXP(-AN7«RB) 

DO 923 KArl ,6 
J-KA+3 

DTEMP(J)zOMEGE(KA)-TEMP(J) 

DTEMP(J+6)-CE(KA)-TEMP(J+6) 

923 CONTINUE 

IF(IPR.EQ.O) 60 TO 12 
IF(NELPP.NE .N) GO TO 12 
IF(NN.NE.NPRIN) GO TO 12 
IF(NCYCLE.EQ.O) -NWALK^D 
NwALK~N WALK+ 1 
NQ-NWALK-Z^J'NCYCLE 
NQQ=NCYCLE -1 
' wRITE{6,2Q) INC 

20 FORMAT (’• INCREMENT’ ,15) 

WRITL(6,750) NIT 
75G FORMAT(’ I TE R A T I ON S ’ , I 5 ) 
wRITE(6,753)N,NN 

753 FORMAT(’ E LE ME N T ’ , I 5 , ' INTEGRATION POINT*, 15) 

IF(Nw.EO.D) i^^RITE(6,23) NQQ 
IF(Ng.GT.O) WRITE(6,39) NCYCLE 

23 FORMAT( * VALUES OF PARAMETERS DURING SOLUTION OF RECYCLE NUMBER* 

1 , 15) 

39 FORMAT(* VALUES OF PARAMETERS DURING ASSEMBLY OF RECYCLE NUMBER* 

1 , 15) 

WRITE (6 •751 ) TEST 

751 FORMAT( ^ ITERATION TEST*,F1C,6> 

WRITE{6,29) 

29 FORMAK* STRAIN INCREMENTS*) 

WRITE(6,30) (DE(J),J-1,NGENS) 

3G FORMAT( 1P6E15.6) 

WRITE(6,31) 

31 FORMAT(* STRESS INCREMENTS*) 

WRITE (6 ,30) (DS ( J ) • J-1 .NGENS ) 

12 RETURN 
END 
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